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SUMMARY 

Central  peaks  are  common  features  observed  in  craters  on  the  earth,  the 
moon,  Mars,  and  Mercury.  Since  these  peaks  do  not  occur  in  all  craters,  they 
should  be  useful  in  providing  strong  constraints  on  both  planetary  evolution 
and  numerical  cratering  simulations.  Unfortunately,  because  the  mechanics  of 
central  peak  formation  has  been  poorly  understood,  little  use  of  those  con- 
straints has  yet  been  made. 

Therefore,  a program  of  numerical  simulations  of  the  ground  response  to  a 
high-explosive  detonation  was  accomplished  to  examine  the  influence  of  model 
conditions  on  calculated  central  peak  formation.  During  this  program,  data  from 
a numerical  simulation  of  a high-explosive  detonation  were  used  as  a surface- 
boundary condition,  and  the  ground  response  was  simulated  by  a computer  code 
that  modeled  two-dimensional,  axisymmetric  problems  of  continuum-mechanics  with 
elastic-plastic  material  models.  First,  a calculation  that  modeled  the  20-ton 
high-explosive  detonation  designated  MIXED  COMPANY  II  showed  that,  when  ball is- 
tically  extrapolated,  the  computed  motions  at  a simulated  time  of  16.4  msec  were 
consistent  with  the  observed  crater  and  formation  of  a central  mound.  The  re- 
sults of  a series  of  calculations  in  which  compaction,  layering,  and  material- 
yield  models  were  varied  indicated  (!)  the  calculated  upward  motions  below  the 
crater  were  eliminated  by  increased  material  compactibi 1 ity,  (2)  the  model  of 
test-site  layering  in  the  MIXED  COMPANY  II  numerical  simulation  only  slightly 
influenced  the  upward  velocities  below  the  crater,  (3)  plastic  volumetric 
increases  of  material  during  Mohr-Coulomb  yield  contributed  significantly  to 
the  calculation  of  upward  motions,  (4)  upward  velocities  for  points  on  the  axis 
of  symmetry  were  first  calculated  where  strength  effects  were  important,  and 
(5)  the  inclusion  of  a lower,  "fluid"  layer  modified  the  calculated  response  in 
an  overlying,  solid  layer  in  a manner  that  may  have  eventually  resulted  in  up- 
ward motions. 

A mechanical  model  of  central  mound  formation  was  developed  with  the  results 
of  the  numerical  calculations  as  a guide.  Material  rebound  in  the  region  where 
strength  effects  are  important  was  emphasized  in  the  model.  Central  mounds 
would  be  inhibited  by  material  compaction  unless  a lower  layer  responded  as  a 
fluid.  The  mechanical  model  includes  enhancement  of  the  central  mound  primarily 
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by  material  bulking  but  also  by  reflections  of  stress  waves  and  the  effect  of 
the  main  shear  wave.  Gravitational  adjustments  that  contribute  to  inward  dis- 
placements are  considered  possible.  This  model  is  found  to  be  consistent  with 
both  the  observed  occurrence  and  structural  relations  of  central  peaks  at  sites 
of  nuclear  and  high-explosive  detonations  and  hypervelocity  impact  events.  The 
conclusions  are  that  the  mechanical  model  is  generally  applicable  to  central 
peak  formation,  the  occurrence  of  a central  peak  in  a crater  is  primarily 
dependent  on  material  properties  of  the  medium,  and  the  calculational  code  used 
for  the  numerical  simulations  can  serve  as  a tool  to  investigate  those  material 
properties . 
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SECTION  I 
INTRODUCTION 

Central  peaks,  or  mounds,  are  a common  feature  of  craters.  Such  peaks  have 
been  observed  in  craters  measured  in  feet  (Roddy,  1968:  1973)  and  in  craters 
measured  in  tens  of  miles  (Baldwin,  1963).  They  occur  in  craters  produced  by 
chemical  explosives  (Roddy,  1973)  and  in  ancient  impact  structures  on  the  earth 
(Howard  et  al . , 1972;  Roddy,  1968;  Dence,  1968;  Beals,  1965).  They  have  been 
seen  in  craters  on  the  Moon  (Baldwin,  1963),  on  Mars  (Hartmann,  1973),  and  on 
Mercury  (Murray  et  al . , 1974).  However,  while  several  authors  have  advanced 
hypotheses  as  to  the  cause  of  central  peaks  (Baldwin,  1963;  Short,  1965; 

Dence,  1968;  Milton  and  Roddy,  1972),  a satisfactory  explanation  of  the  mechan- 
ics of  central  mound  formation  has  not  been  demonstrated. 

Central  mounds  serve  as  a very  useful  constraint  on  cratering  calculations 
because  they  are  directly  observable.  The  Defense  Nuclear  Agency  (a  branch  of 
the  U.S.  Department  of  Defense)  has  sponsored  many  attempts  to  simulate  numeri- 
cally the  cratering  and  ground-shock  effects  of  experiments  where  high-explosive 
charges  were  detonated  at  the  surface  of  the  earth  (Christensen,  1969;  Maxwell 
and  Mcises,  1971b;  Port  and  Gajewski,  1973;  Wright  et  al.,  1973;  Ialongo,  1973). 
Primary  constraints  on  these  calculations  have  been  comparisons  with  crater 
radii  and  depths,  and  comparisons  with  ground-motion  data  obtained  from  active 
instrumentation  usually  emplaced  more  than  one  crater  radius  from  ground  zero. 
Although  central  mounds  occurred  in  many  of  the  craters  (Roddy,  1968;  1973), 
little  attention  has  been  given  to  whether  the  motions  predicted  by  such  calcu- 
lations are  consistent  with  central  peak  formation.  Christensen  (1970)  con- 
sidered such  a constraint,  and  he  found  the  requirements  for  upward  motions 
below  the  detonation  point  to  be  a very  useful  guide  to  the  physics  required 
in  the  calculation  because,  while  all  computations  predicted  a crater,  a cen- 
tral mound  was  predicted  only  under  limited  conditions.  Therefore,  the  presence 
or  absence  of  a central  mound  in  an  experimentally  produced  crater  provides  an 
immediate  direct  test  of  the  adequacy  of  numerical  simulations,  and  an  under- 
standing of  the  causes  of  the  central  mound  will  be  important  to  the  under- 
standing of  cratering  and  ground-shock  mechanics. 
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Also,  the  centra]  peaks  observed  in  hypervelocity  impact  craters  on 
planetary  surfaces  provide  information  on  the  conditions  of  the  impact  event. 

The  high-velocity  impact  of  a body  on  a planetary  surface  is  a dynamic  test  of 
the  two  bodies  over  a very  short  time  period.  This  information  is  recorded  in 
the  occurrence  of  a central  mound;  because,  while  the  size  of  the  crater  is 
primarily  controlled  by  the  mass  and  velocity  of  the  impacting  particle 
(Baldwin,  1963),  the  formation  of  a central  mound  is  primarily  controlled  either 
by  properties  of  the  planet  (Baldwin,  1963;  Short,  1965;  Dence,  1968),  by 
properties  of  the  impacting  body  (Milton  and  Roddy,  1972),  or  by  properties  of 
both  (Roddy,  1968b).  An  understanding  of  the  causes  of  central  mound  formation 
may  allow  us  to  obtain  information  concerning  the  conditions  at  the  time  a 
crater  was  formed.  The  evolution  of  those  conditions  may  then  be  studied  by 
observations  of  craters  of  different  ages,  and  variations  of  those  conditions 
with  location  may  be  tested  with  central-peak  data  in  different  locations. 

The  central  peaks  of  craters  represent  important  sampling  sites  for  any 
extraterrestrial  landing  or  remote  sensing  mission.  Roddy  (1968)  showed  that 
the  material  in  the  central  mound  is  the  deepest  material  exposed  during  the 
cratering  event.  A traverse  across  the  mound  will  sample  the  deepest  strati- 
graphic section  obtainable  at  the  surface  of  a crater.  An  understanding  of 
central  peak  formation  will  aid  the  determination  of  the  pre-impact  location 
of  the  material . 

Therefore,  a program  of  numerical  simulations  was  accomplished  to  examine 
the  causes  of  central  peak  formation  in  shock-wave  cratering  event.  During  this 
program,  the  models  were  limited  to  the  simplest  possible  expressions  to  demon- 
strate which  factors  were  most  important  in  the  formation  of  a central  mound. 
These  models  were  an  idealization  of  much  more  complex  material  behavior.  The 
AFT0N-2A  (see  appendix  B)  computer  code  was  used  because  it  was  already  active 
in  ground-shock  calculations  (see,  for  example,  Port  and  Gajewski,  1973)  and 
because  no  code  development  was  required  other  than  modifications  of  material 
behavior  models.  One  model  of  a high-explosive  detonation  was  used  for  all  of 
the  numerical  simulations,  and  only  models  of  material  properties  were  varied. 
The  results  of  the  calculations  were  generalized  to  different  types  of 
cratering  events  by  developing  a general  model  for  the  causes  of  central  mound 
formation  that  was  compared  to  reported  observations  of  central  mounds. 

To  clarify  this  discussion,  the  terms  "shock-wave  cratering  event"  and 
"central  mount"  should  be  defined.  A "shock-wave  cratering  event"  is  an  event 
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that  transfers  a large  amount  of  energy  to  a small  volume  of  a halfspace  by 
sending  a shock  wave  into  the  halfspace  and  forms  a crater  at  the  surface.  This 
term  applies  to  both  a hypervelocity  impact  and  a surface  chemical  or  nuclear 
explosive  detonation.  Further,  as  Shoemaker  (1961)  suggests,  this  term  is 
more  basic  than  explosion  cratering  event  because,  even  if  the  characteristics 
of  craters  are  controlled  by  the  expansion  of  gases  near  the  source  of  the 
event  (Baldwin,  1963),  the  passage  of  the  shock  wave  through  the  material  is 
the  mechanism  which  establishes  the  conditions  for  such  expansion.  A "central 
mound  is  a local  topographical  high  at  the  center  of  a crater  that  is  composed 
of  material  that  was  displaced  upward  during  the  cratering  event.  This  term 
is  applicable  to  a definite  structural  feature,  and  is  not  meant  to  include 
the  possibility  that  material  ejected  from  the  crater  may  subsequently  fall 
into  the  crater  and  form  a hill  at  the  center.  Also,  in  this  report,  the 
terms  "central  peak"  and  "central  mound"  will  be  synonymous. 

Three  assumptions  were  basic  to  the  physical  models  used.  First,  all  calcu- 
lations were  performed  assuming  axial  symmetry,  i.e.,  all  properties  can  be 
described  in  terms  of  the  radial  and  axial  position  coordinates  of  a cylindrical 
coordinate  system.  Second,  all  material  models  were  assumed  to  be  isotropic. 
Finally,  no  energy  transfer  by  radiation  or  conduction  was  considered. 

In  addition  to  the  basic  assumptions,  several  additional  assumptions  were 
involved  in  this  study.  Some  of  these  assumptions  were  inherent  in  the  AFT0N-2A 
code  and  are  described  in  appendix  B.  Other  assumptions  were  involved  in  the 
description  of  the  numerical  experiments  accomplished  and  are  described  during 
the  presentation  of  those  problems. 

The  generalization  of  the  results  of  numerical  simulations  should  be  done 
cautiously  because  those  results  strictly  apply  only  to  definite  models  with 
specific  input  parameters.  The  following  procedure  was  used  to  generalize 
specific  numerical  results  concerning  the  mechanics  of  central  peak  formation. 
First,  the  available  information  pertinent  to  central  peaks  was  reviewed  to 
provide  the  broadest  possible  base  of  data.  Then,  a numerical  model  of  the 
high-explosive  experiment  MIXED  COMPANY  II  is  described  to  demonstrate  the 
applicability  of  the  numerical  results  to  that  one  experiment.  The  results 
from  additional  numerical  simulations  in  which  material  models  were  varied  are 
described  to  determine  what  properties  are  important  to  central  peak  formation. 
These  properties  are  included  in  a general  model  of  central  peak  mechanics  which 
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SECTION  II 

PREVIOUS  WORK  OF  OTHERS 

Three  types  of  information  are  available  concerning  central  mounds  in 
shock-produced  craters.  First,  observations,  presented  by  others,  of  the 
occurrence  of  central  peaks  and  of  the  structure  of  central  uplifts  in  both 
ancient  impact  features  on  earth  and  experimental  high-explosive  detonations 
provide  constraints  on  any  explanation  of  central  peak  formation.  Second, 
numerical  simulations  of  shock-wave  cratering  events  provide  guides  to  the 
physical  processes  involved  in  central  mound  formation.  Finally,  several  others 
have  previously  proposed  hypotheses  concerning  the  causes  of  central  peak  for- 
mation. In  this  section,  information  of  each  type  is  reviewed  to  provide  a 
basis  for  later  conclusions  concerning  the  mechanics  of  central  mound  formation. 

1 . OCCURRENCE  AND  STRUCTURE  OF  CENTRAL  PEAKS 

Observations  of  the  structure  of  central  uplifts  and  their  occurrence  in 
craters  indicate  that  similar  relations  apply  to  both  hypervelocity  impact 
events  and  explosive  detonations.  The  material  in  the  central  uplifts  of  both 
ancient  impact  structures  on  earth  and  high-explosive  craters  is  displaced 
upward  from  its  original  position  (Roddy,  1968).  Horizontal  displacements  of 
the  material  that  forms  central  mounds  are  probably  inward  in  the  deeper  regions 
and  outward  in  the  shallower  regions  (Howard,  et  al.,  1972;  Milton,  et  al., 
1972).  Shatter  cones  are  frequently  found  in  the  central  uplifts  (Dietz  1968; 
Roddy,  1973),  indicating  that  maximum  stresses  were  on  the  order  of,  but 
above,  the  Hugoniot  elastic  limit  of  the  material.  Central  mounds  do  not  occur 
in  craters  formed  in  very  porous  material. 

a.  Central  Peaks  in  Hyperve iocity  Impact  Structures 

Central  uplifts  have  been  observed  in  many  structures  that  have  effects 
which  are  commonly  associated  with  sites  of  hypervelocity  impact  events  on 
earth  (Dietz,  1968).  Roddy  (1968)  examined  one  of  those  structures  (2.3  miles 
in  diameter)  at  Flynn  Creek,  Tennessee  (figure  1).  He  found  that  Stones  River 
and  Knox  strata  occur  in  the  center  as  folded,  faulted,  and  brecciated  material 
which  forms  a hill  with  a top  about  370  feet  above  the  original  crater  floor. 
This  hill  had  the  general  structure  of  a domed  megabreccia  block  with  100  feet 
of  the  Knox  formation,  now  exposed  as  steeply  dipping  strata,  raised  1100  feet 
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above  its  original  position.  Where  the  base  of  the  mound  intersected  the 
bedded  breccia  on  the  crater  floor,  the  hill  was  almost  3000  feet  in  diameter 
with  sides  sloping  an  average  of  15°.  Similar  upward  displacements  of  the 
material  in  central  mounds  were  observed  at  the  Wells  Creek  structure, 

Tennessee  (Stearns,  et  al . , 1968);  Sierra  Madera,  Texas  (Howard,  et  al.,  1972); 
and  Gosses  Bluff,  Australia  (Milton,  et  al . , 1972). 

Some  evidence  of  inward  displacements  at  depths  also  existed  in  many 
of  the  structures.  Howard,  et  al.,  (1972)  suggested  that  individual  beds  in 
the  Sierra  Madera  uplift  were  faulted  and  folded  to  an  extent  that  the  total 
strike  length  of  each  bed  was  greater  than  the  perimeter  on  which  it  lies. 

This  shortening  may  have  been  as  great  as  25  percent  in  some  stratigraphic 
sections,  although  this  estimate  was  based  on  the  possibly  invalid  (Milton  and 
Roddy,  1972)  assumption  that  the  displacement  of  each  segment  could  be  resolved 
into  translation  plus  rotation  about  no  more  than  one  axis.  Also,  stratigraphic 
beds  appeared  to  have  been  thickened  such  that  near-vertical  beds,  which  were 
no  more  than  1200  feet  thick,  filled  a minimum  width  of  5000  feet. 

The  orientation  of  shatter  cones  (shatter  cones  will  be  discussed  later) 
has  also  been  used  as  a measure  of  inward  displacement  (Milton,  et  al . , 1972). 
Measurements  of  such  orientations  at  Gosses  Bluff  indicated  that  when  a common 
shock  focus  was  assumed  inward  displacements  were  from  20  to  52  percent  of  the 
original  radial  distance  from  the  center,  with  the  deeper  strata  displaced 
inward  more  than  the  shallower  strata.  However,  these  estimates  could  be 
significantly  reduced  if  the  shock  was  produced  by  a vertical  line  source. 

A complete  elimination  of  inward  displacements  would  require  that  the  line 
source  was  6300  feet  long.  If  the  relation  that  upper  heds  moved  inward  rela- 
tive to  lower,  as  suggested  by  Howard,  et  al . , (1972)  based  on  fold  patterns 
at  Sierra  Madera,  is  also  valid  at  Gosses  Bluff,  then  the  assumption  of  a 
common  source  results  in  an  incorrect  relation  of  displacements  between  strata, 
and  quantitative  estimates  that  are  based  on  that  assumption  are  not  valid. 

While  evidence  of  inward  displacements  in  the  deeper  regions  is  not 
complete,  outward  displacements  in  much  shallower  regions  are  observed.  For 
example,  at  Gosses  Bluff  the  upper  ends  of  layers  lie  as  overturned  plates  or 
detached  blocks  on  the  truncated  edge  of  stratigraphical ly  higher  units 
(Milton,  et  al . , 1972).  In  addition,  330  foot  long  blocks  of  sandstone  lie 
1000  feet  from  their  stratigraphic  outcrop,  indicating  an  outward  ballistic 
flight. 
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The  peak  shock  pressure?  experienced  by  the  material  in  central  mounds 
can  be  estimated  on  the  basis  of  the  occurrence  of  shock  effects.  One  macro- 
scopic shock  effect,  already  mentioned,  is  the  shatter  cone.  Dietz  (1968) 
describes  shatter  cones  as  cup-and-cone  structures  with  striated  surfaces 
that  radiate  from  small  half-cones  on  the  face  of  a master  cone.  They  are 
most  common  in  carbonate  rocks,  but  are  also  known  in  shale,  sandstone,  quartz- 
ite, and  other  lithologies.  They  frequently  occur  in  central  uplifts,  as 
shown  in  table  1,  indicating  that  conditions  favorable  to  shatter-cone  formation 
are  experienced  in  that  material.  A theoretical  study  of  shatter-cc.iing  (as 
reported  by  Dietz,  1968)  shows  that  shatter  cones  are  shock  fractures  formed 
along  a travelling  boundary  between  the  plastic  and  elastic  response  of  a 
material  defined  by  the  dynamic  elastic  limit,  with  the  plastic  domain  moving 
relative  to  an  elastic  domain.  The  analysis  is  consistent  with  the  observations 
that  shatter  cones  appear  to  be  formed  prior  to  significant  material  displace- 
ment, and  high-pressure  phases  (coesite,  stishovite,  maskelynite)  have  not  been 
found  associated  with  shatter  cones.  Thus,  shatter  cones,  and  by  association 
central  mounds,  appear  to  be  formed  in  material  where  the  shock  pressures  were 
close  to,  but  above,  the  Hugoniot  elastic  limit. 


Table  1 

STRUCTURES  WITH  SHATTER  CONES  IN  CENTRAL  UPLIFT 
(Dietz,  1968) 


Structure 

Location 

Rock  type  of  shatter  cones 

1.  Steinheim  Basin 

Germany 

Limestones 

2.  Wells  Creek  Basin 

Tennessee 

Dolomite 

3.  Crooked  Creek 

Mi  ssour* 

Limestone 

4.  Serpent  Mound 

Ohio 

Limestone 

5.  Flynn  Creek 

Tennessee 

Limestone 

6.  Sierra  Madera 

Texas 

Limestone 

7.  Verdefort  Ring 

South  Africa 

Granite 

8.  Clearwater  Lake  West 

Quebec,  Canada 

Granitic  Gneiss 

9.  Sudbury 

Ontario,  Canada 

Quartzite,  Shale,  Granite 

10.  Manicouagan-Mushalagan 

Quebec,  Canada 

Crystalline  Gneisses 

1 1 . Gosses  Bluff 

Austral ia 

Limestone,  Sandstone,  etc. 
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An  extensive  study  of  circular  structures  that  have  effects  consistent 
with  hypervelocity  impact  sites  has  led  Dence,  et  al.,  (1968)  to  conclude 
that  there  is  a critical  crater  size  required  to  form  central  uplifts  in 
granitic  gneisses.  Information  from  structural  mapping,  gravity  surveys, 
and  drilling,  combined  with  evidence  of  previous  shock-wave  experience  in 
material s , has  indicated  that  there  are  at  least  12  shock-produced  craters  on 
the  Canadian  Shield.  The  smaller  of  these  craters,  with  diameters  of  2.5  miles 
or  less  (Brent,  Holleford,  New  Quebec  craters),  have  a bowl-shaped  structure 
with  no  central  uplift.  The  craters  with  diameters  greater  than  5.5  miles 
(Deep  Bay,  Clearwater  Lakes,  Carswell  Lake),  however,  show  a complex  structure 
which  includes  a central  uplift,  an  annulus  of  brecciated  rock,  arid  a peripheral 
depression  which  surrounds  the  crater. 

Observations  of  the  occurrence  of  central  mounds  on  the  Earth,  Mars, 
and  the  Moon  have  been  interpreted  to  show  that  gravity  has  an  influence  on 
the  occurrence  of  central  peaks  (Hartmann  1972,  1973).  Hartmann  (1972)  sug- 
gested that  data  on  the  size  distribution  of  craters  with  central  peaks  as  a 
function  of  crater  diameter  (figure  2)  indicated  that  central  peaks  tended  to 


Figure  2.  The  Size  Distribution  of  Craters  with  Central  Peaks  on  the 
Earth,  Mars,  and  the  Moon  (Solid  line  on  earth  data  includes 
structural  uplifts  in  astroblemes  and  is  based  on  a total  of 
33  structures  (after  Hartmann,  1972) 
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form  in  craters  of  smaller  diameter  as  gravity  increased;  although  the  statisti- 
cal base  for  the  Earth  data  was  only  33  cryptovolcanic  structures.  He  suggested 
the  relationship  was 


D oc  q“l*25t0.2 

c 3 

where  g was  the  gravitational  acceleration  and  Dc  represented  either  the  minimum 
diameter  for  craters  with  central  peaks  or  the  diameter  of  craters  with  the 
maximum  frequency  of  central  peak  occurrence.  The  data  also  showed,  however, 
that  the  inferred  minimum  diameter  of  craters  with  central  peaks  ^as  signifi- 
cantly different  from  the  diameter  at  maximum  frequency.  This  difference 
indicated  that  gravitational  stress,  while  a contributory  factor,  was  not  the 
only  cause  of  central  mound  formation. 

b.  Central  Mounds  in  High  Explosive  Detonation  Experiments 

Similar  structural  relations  have  been  observed  at  central  mounds  in 
craters  caused  by  large  chemical-explosive  detonations.  One  series  of  tests 
was  located  at  the  Watching  Hill  Test  Range  near  Suffield,  Alberta,  and  was 
sponsored  by  the  Defense  Research  Establishment  of  Canada.  Programs  during 
this  series  included  SNOWBALL,  DISTANT  PLAIN,  PRAIRIE  FLAT,  and  DIAL  PACK. 

The  geology  at  the  test  range  was  characterized  by  a ground  water  table  usually 
near  25  feet  depth  (Zelasko  and  Baladi,  1971).  The  presence  and  depth  of  this 
water  table  resulted  in  essentially  a two-layer  structure.  The  material  above 
25  feet  depth  was  a low-density  soil  that  displayed  a Mohr-Coulomb  yield  surface 
with  a slope  near  1.  The  material  below  25  feet  depth  was  a denser,  saturated 
soil  that,  for  confining  pressures  less  than  40  bars,  had  little  strength 
dependence  on  confining  pressure.  Also,  the  Poisson's  ratio  of  the  material 
increased  from  0.30  near  the  surface  to  0.47  at  30  feet  depth. 

Central  peaks  with  characteristics  similar  to  the  central  uplifts 
observed  in  ancient  impact  structures  were  formed  in  many  of  the  craters  that 
resulted  from  the  explosive  tests.  Roddy  (1968)  described  the  crater  (figure  3) 
from  the  500-ton  TNT  event,  SNOWBALL,  as  shallow  and  flat-floored  with  a diam- 
eter of  more  than  300  feet  and  a maximum  depth  of  22  feet.  The  central  mound, 
which  was  nearly  19  feet  high,  consisted  of  folded  and  faulted  clay  beds  in  a 
tightly  folded  dome.  The  beds  showed  plastic  thickening  and  thinning  with  a 
stratigraphic  horizon  lifted  nearly  24  feet.  During  the  test  series  at  Suffield, 
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Figure  3.  Cross  Section  of  the  Crater  Produced  by  the  500-Ton 
TNT  Event  SNOWBALL  (after  Roddy,  1968) 

a 20-ton  spherical  charge  and  a 100-ton  hemispherical  charge  formed  craters  with 
no  central  peaks,  while  a 100-ton  spherical  charge  and  500-ton  hemispherical  and 
spherical  charges  formed  craters  with  central  mounds  (Milton  and  Roddy,  1972). 

A series  of  high-explosive  detonations,  named  MIDDLE  GUST,  was  performed 
near  Pueblo,  Colorado,  during  1971  and  1972  (Myers,  1973).  This  series  included 
a total  of  five  experiments  at  two  sites  (see  table  2).  One  site  had  a 10-foot 

Table  2 


DEFENSE  NUCLEAR  AGENCY  HIGH  EXPLOSIVE  EXPERIMENTS 


Experiment 

Charge 

Charge 

Central 

(tons  TNT) 

position 

Test  side 

mound 

MIDDLE  GUST  I 

20 

hal  f-buried 

"wet"  shale 

2 ft 

MIDDLE  GUST  II 

100 

elevated 

"wet"  shale 

5 ft 

MIDDLE  GUST  III 

100 

surface 

tangent 

"wet"  shale 

trough 

MIDDLE  GUST  IV 

100 

surface 

tangent 

"dry"  shale 

2 ft 

MIDDLE  GUST  V 

20 

hal f-buried 

"dry"  shale 

3 ft 

MIXED  COMPANY  I 

20 

half-buried 

sandstone 

3 ft 

MIXED  COMPANY  II 

20 

surface 

tangent 

sandstone 

7 ft 

MIXED  COMPANY  III 

500 

surface 

tangent 

sandstone 

5 ft 
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overburden  of  sandy  clay  over  fractured  clay  shale  that  interfaced  with  com- 
petent shale  23  feet  below  the  ground  surface  (Windham,  et  al . , 1973).  This 
site  was  called  the  "wet"  site  because  a perched  water  table  on  the  top  of 
the  :ompetent  shale  extended  to  about  4 feet  below  the  ground  surface.  Two 
sets  of  nearly  vertical  joints  existed  in  the  competent  shale.  The  sets  were 
nearly  perpendicular  and  had  intervals  between  joints  of  6 to  8 feet  and  10 
to  14  feet.  The  second  or  "dry"  site  was  23  feet  of  fractured  clay  shale 
over  a more  competent  shale  with  no  near-surface  water  table.  Three  sets  of 
nearly  vertical  joints  existed  in  the  weathered  shale  at  the  second  site.  The 
maximum  strength  of  even  the  competent  shale  at  both  sites  was  less  than  100 
bars  and  independent  of  confining  pressure.  All  the  craters  extended  into  the 
weathered  shale,  and  the  MIDDLE  GUST  III  crater  reached  the  competent  shale 
(Myers,  1973).  All  of  the  craters  except  MIDDLE  GUST  III  had  interior  mounds 
that  were  2 to  5 feet  high,  although  the  top  of  one  of  the  mounds  was  offset 
from  ground  zero.  The  mounds  tended  to  fracture  along  old  joint  directions 
(Roddy,  1973).  The  MIDDLE  GUST  III  crater  had  a 5-foot  deep  central  trough 
in  the  competent  shale.  Roddy  (1973)  reported  shatter  cones  in  the  MIDDLE  GUST 
IV  and  V central  mounds. 

An  additional  series  of  three  high  explosive  detonations,  called  the 
MIXED  COMPANY  series,  was  performed  near  Grand  Junction,  Colorado,  during 
1972  (Choromokos  and  Kelso,  1973).  The  sites  for  these  experiments  were  surface 
layers  of  alluvial  sandy  soil  over  sandstone  with  no  significant  water  content. 
The  alluvial  soil  layer  was  5 feet  thick  for  events  I and  III  and  1.8  feet 
thick  (Day,  1973)  for  event  II.  The  sandstone  was  generally  weathered  to  a 
depth  of  12  feet  below  the  surface. 

The  craters  that  resulted  from  these  experiments  all  had  central  uplifts 
(Roddy,  1973).  The  first  event  produced  a crater  with  an  apparent  depth  of 
15  feet  and  a central  peak  3 feet  high.  The  sandstone  beds  in  this  uplift  were 
generally  intact  on  the  flanks  but  were  brecciated  in  the  core.  The  second 
event  produced  a crater  with  an  apparent  depth  of  7 feet.  A very  large  central 
mound  covered  the  crater  floor  and  extended  nearly  3 feet  above  the  original 
ground  level.  The  third  event  produced  a crater  with  an  apparent  depth  of 
18  feet  and  a poorly  formed  irregular  dome  of  massive  sandstone  5 feet  high. 

The  crater  floor  surrounding  the  uplift  consisted  of  large  slabs  of  sandstone 
that  sloped  upward  towards  ground  zero  and  exhibited  both  fracturing  and 
faulting  approximately  parallel  to  the  local  joint  pattern.  A circular  ring 
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fault  formed  on  the  crater  floor  at  the  base  of  the  crater  walls,  locally 
separating  the  floor  from  the  walls.  Carnes  (1973a)  reported  that  permanent 
displacements  in  the  sandstone  beyond  the  crater  walls  were  primarily  upward 
and  outward. 

c.  Craters  Without  Central  Peaks 

Central  mounds  were  not  observed  in  all  shock-produced  craters.  In 
particular,  there  was  a notable  lack  of  central  peaks  produced  by  nuclear 
detonations  at  the  Nevada  Test  Site  (Roddy,  1968)  and  at  Eniwetok  and  Bikini 
Atolls  in  the  Pacific  (Circeo  and  Nordyke,  1964).  Such  nuclear  experiments 
included  detonations  with  yields  from  a kiloton  to  over  10  megatons  and 
deeply  buried,  near  surface,  and  above  surface  shotpoint  locations.  The  test 
sites  were  dry  alluvium  at  the  Nevada  Test  Site  and  unconsolidated  sands 
and  gravels  over  coral  reefs  in  the  Pacific.  The  common  geologic  characteris- 
tic of  these  sites  is  the  porous  structure  of  the  cratered  materials.  The 
medium  beneath  the  Barringer  Crater  in  Arizona,  a meteorite  impact  crater 
with  no  central  mound,  is  porous  sandstone  (Shoemaker,  1963).  Also,  craters 
without  central  mounds  were  produced  in  impact  craters  formed  in  dry,  non- 
cohesive  quartz  sand  during  laboratory  tests  (Gault,  et  al . , 1968). 

2.  PREVIOUS  CALCULATIONAL  PROGRAMS 

Several  attempts  have  been  made  by  others  to  numerically  simulate  shock- 
wave  cratering  events.  These  attempts  have  shown  that  the  ground-motion  history 
during  a simulation  was  sensitive  to  the  amount  of  material  compaction  that 
was  modeled  for  a complete  cycle  of  stress  loading  and  unloading.  A second 
important  feature  of  the  simulations  was  that  the  motion  caused  by  the  simulated 
shear  wave  was  toward  the  axis  of  symmetry. 

a.  Distant  Plain  6 Simulation 

The  coupled  Eulerian-Lagrangian  computer  code  called  ELK  was  used  for 
three  attempts  to  numerically  simulate  the  crater  and  central  mound  produced 
by  the  100-ton,  surface-tangent  high-explosive  experiment.  Distant  Plain  6 
(Christensen,  et  al.,  1968;  Christensen,  1970).  The  material  models  for  all 
three  of  the  calculations  were  based  on  reported  test-site  data;  however,  only 
one  of  the  simulations,  ELK  31,  included  a precompaction  model  beneath,  and 
as  a result  of,  the  100-ton  explosive  charge.  The  ELK  31  calculation  was 
continued  until  a simulated  time  of  220  msec  (figure  4)  and  showed  the  develop- 
ment of  upward  velocities  near  the  axis  of  symmetry  after  160  msec.  At  220 
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msec,  the  flow  pattern  showed  that  material  near  the  crater  boundary  at  the 
20-foot  range  was  moving  down  and  toward  the  vertical  axis.  This  motion  re- 
sulted in  a vortex  pattern  centered  at  a range  of  12  feet.  Extension  of  this 
flow  pattern,  accomplished  by  extrapolating  deceleration,  resulted  in  calculated 
crater  dimensions  at  a simulated  time  of  1 second  that  were  consistent  with 
the  observed  crater.  The  other  two  calculations  were  stopped  by  a simulated 
time  of  125  msec  because  the  computed  depths  of  the  craters  were  too  great. 
Christensen  (1970)  concluded  that  causes  for  the  upthrust  included  the  airblast- 
induced  shear  wave  interacting  at  the  axis  of  symmetry,  the  effect  of  gravity, 
and  the  compaction  cone  that  resulted  from  quick  settlement  under  the  100-ton 
explosive  load,  but  he  did  not  determine  the  contribution  of  each  cause. 

b.  Mine  Under  Simulation 

A numerical  simulation  was  also  accomplished  to  model  the  test  event 
MINE  UNDER  (Maxwell  and  Moises,  1971a).  The  event  was  the  detonation  of  a 
100-ton  spherical  charge  of  TNT  over  granite.  The  charge  was  centered  at  a 
height  of  2 charge  radii  to  produce  only  airblast  loading  on  the  ground.  The 
test  site  was  composed  of  weathered  granite  with  a compressional  wave  speed  of 
10,000  ft/sec  and  a shear  wave  speed  of  6100  ft/sec.  Sample  porosities  varied 
a factor  of  2 from  a mean  value  of  5 percent.  This  porosity  resulted  in  a 
residual  compression,  after  a cycle  of  compressive  loading  and  unloading,  of 
20  percent  of  the  peak  compression  for  peak  pressures  below  43  kbar.  A complex 
yield  model  which  included  brittle  fracture  and  sliding  on  cracks  was  also 
included  in  the  calculation,  with  a von  Mises  limit  of  30  kbar  reached  by  a 
pressure  of  24  kbar.  The  results  of  this  calculation  showed  that  the  calculated 
shear  wave,  supported  by  the  strong  rock  model,  caused  a clockwise  rotation  in 
the  material  flow  pattern.  However,  data  from  instrumentation  in  the  actual 
event  did  not  indicate  such  a substantial  shear  wave.  Maxwell  and  Moises 
concluded  that  the  in  situ  rock  strength  was  much  lower  than  the  strength 
included  in  the  model. 

c.  Sierra  Madera  Simulation 

A numerical  simulation  of  the  event  which  may  have  formed  the  Sierra 
Madera  formation  was  also  accomplished  by  Maxwell  and  Moises  (1971b).  For  this 
simulation  a sphere  with  a radius  of  328  feet  and  a velocity  of  19  miles/sec 
was  assumed  to  impact  vertically  on  a halfspace.  Both  the  sphere  and  the  half- 
space were  assumed  to  be  composed  of  the  same  material,  which  had  an  assumed 
density  of  2.7  gm/cc.  The  parameters  of  the  .material  equation  of  state  were 
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based  only  on  Hugoniot  data  for  basalt.  The  yield  model  was  a 0.2-kbar  von 
Mises  limit  until  a calcul ational  zone  experienced  zero  pressure,  after  which 
that  zone  was  assumed  to  have  no  shear  strength.  The  calculation  resulted  in 
upward  velocities  below  the  impact  point  by  a simulated  time  of  5.5  seconds 
with  a toroidal  flow  pattern  developed  by  9.5  seconds  that  continued  until 
the  calculation  was  terminated  at  30  seconds.  Maxwell  and  Moises  concluded 
that  the  dominant  driving  force  of  the  central  uplift  was  the  release  of  the 
overburden  by  excavation.  The  entire  flow  pattern  after  5.5  seconds,  however, 
could  be  explained  by  the  flow  of  a liquid  under  the  influence  of  gravity. 

d.  Nuclear  Explosion  Simulation 

A series  of  calculations,  called  REVROC,  was  completed  at  the  Air  Force 
Weapons  Laboratory  (AFWL,  1973)  to  study  the  effects  of  layered  bedrock  on 
calculated  near-surface  ground  motions  caused  by  a simulated  nuclear  explosion. 
A two-dimensional,  axisymmetric  computer  code  was  used  for  the  simulation. 
Material  models  included  irreversible  compaction  after  a cycle  of  loading  and 
unloading.  The  results  showed  that,  by  0.5  second,  the  calculated  flow  field 
included  upward  motion  near  the  vertical  axis  of  symmetry.  These  motions 
seemed  to  be  caused  by  the  primary  shear  wave  and  occurred  even  in  the  bottom 
layer  of  material.  Also,  the  motions  were  a function  of  the  amount  of  irre- 
versible compaction  included  in  the  model,  with  less  compaction  favoring  more 
upward  motion. 

e.  MIDDLE  GUST  III  Simulation 

Port  and  Gajewski  (1973)  performed  a numerical  parametric  study  to 
examine  the  causes  of  major  discrepancies  between  the  ground  motions  calculated 
using  pretest  models  of  the  MIDDLE  GUST  III  event  and  ground  motions  measured 
during  the  experiment.  The  main  discrepancy  was  the  failure  of  the  numerical 
calculations  to  simulate  accurately  the  arrival  and  magnitude  of  the  large 
direct-induced  ground  shock  which  dominated  the  experimental  motion.  Three 
alternate  material  models  were  evaluated  in  their  calculations.  The  first, 
or  laboratory,  model  was  based  on  detailed  laboratory  uniaxial-strain  tests 
on  samples  of  materials  obtained  from  site  drill  cores.  The  second,  or  CIST, 
model  was  based  on  cylindrical  in  situ  test  results  of  the  MIDDLE  GUST  site. 

The  CIST  test  was  used  to  determine  dynamic  moduli  of  in  situ  materials  by 
measuring  the  ground  motions  caused  by  a cylindrically  symmetric  shock  input 
(Davis,  1973).  The  third  model  was  based  on  seismic  velocity  data  of  the  site 
and  was  referred  to  as  the  seismic  model.  The  first  two  models  included 
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irreversible  compaction  after  a load-unload  pressure  cycle,  while  the  seismic 
model  was  incompactible.  The  velocities  of  large-amplitude  stress  waves 
implied  by  the  models  were  lowest  for  the  laboratory  model  and  highest  for 
the  seismic  model.  The  same  equation  of  state  was  used  for  high  pressures  in 
all  three  models. 


The  results  of  the  parametric  study  showed  that  the  laboratory  model 
was  inadequate.  That  mo<Jel  resulted  in  wave  speeds  that  were  one-third  to 
one-fourth  of  the  values  required  to  match  the  arrival  time  of  strong  ground- 
motion  signals.  Further,  that  model  failed  to  produce  the  magnitude  of  the 
peak  upward  velocities  near  the  surface  caused  by  the  direct-induced  wave 
(figure  5).  The  shock  arrival  times  that  were  measured  above  23  ft  depth 
were  most  consistent  with  the  CIST  model  calculation.  The  data  from  instru- 
ments placed  below  23  feet  indicated  wave  speeds  greater  than  even  the  seismic 
model.  The  peak  upward  velocities  near  the  surface  were  matched  only  with  the 
seismic  model.  The  crater  profiles  predicted  for  all  three  models  were  nearly 
the  same,  with  a maximum  depth  below  ground  zero  near  17  feet,  while  the 
observed  depth  was  21  feet. 

3.  POSTULATED  MECHANISMS 

Several  mechanisms  have  been  proposed  to  explain  the  formation  of  central 
peaks.  These  mechanisms  may  be  broadly  divided  into  (1)  effects  related  to 
stress  waves  and  (2)  rapid  gravitationally  controlled  adjustments  of  the  walls 
of  an  initial  crater.  The  first  of  these  broad  divisions  includes  rebound, 
stress  wave  reflections,  shear  wave  effects,  and  special  boundary  conditions 
caused  by  the  impact  of  low-density  bodies.  The  second  division  includes  deep 
gravitational  sliding  and  Rayleigh  jet  formation. 

a.  Stress  Wave  Related  Mechanisms 

The  rebound  of  material  below  the  crater  following  the  compression  by 
the  shock  wave  was  suggested  by  Boon  and  Albritton  (1938)  as  the  primary 
mechanism  for  forming  central  mounds.  This  rebound  results  from  the  accelera- 
tion of  material  toward  the  stress  wave  source  as  a result  of  the  decreasing 
stress  gradient  that  extends  to  the  free  surface.  Baldwin  (1963)  used  a set 
of  two  very  small-scale  explosive  experiments  to  examine  this  mechanism. 

In  the  first,  a 40-grain  dynamite  charge  was  detonated  1 inch  below  the  surface 
of  a specially  built-up  volume  of  soil.  A box  3 feet  square  and  1 foot  deep 
was  filled  with  soil.  The  bottom  6 inches  consisted  of  ordinary  soil.  Above 
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this  soil  were  six  horizontal  layers  of  colored  soil,  each  1-inch  thick,  with 
each  layer  lightly  tamped  into  place.  The  second  test  was  a repeat  of  the 
first  with  the  soil  colored  into  vertical  layers.  From  these  two  tests,  a 
composite  crater  (figure  6)  was  determined.  Baldwin  concluded  that  a flat  area 


Figure  6.  Composite  Crater  Profile  Constructed  from  the 
Detonation  of  Two  40-Grain  Dynamite  Charges 
with  the  Shot  Point  at  the  Base  of  the  Black 
Layer  (The  curved  dotted  line  gives  the  limit 
of  the  volume  from  which  the  soil  was  actually 
blasted  from  the  crater;  the  lens  at  the 
bottom  of  the  crater  is  white  material  of 
lower  than  original  density  (after  Baldwin, 

1963,  pg.  120)). 

in  the  center  of  the  crater  was  an  incipient  central  cone  formed  by  rebound. 
This  conclusion  was  reached  because  some  of  the  horizontal  white  layer  was 
found  above  the  red  layer  and  partially  under  the  gray  layer.  This  white 
material  was  of  lower-than-normal  density,  while  the  yellow  and  white  layers 
below  it  were  denser  than  normal.  Dence  (1968)  also  mentioned  that  unloading 
of  materials  after  the  shock  wave  might  influence  central  mound  formation, 
but  he  did  not  emphasize  this  mechanism  because  of  the  small  increase  in 
specific  volumes  for  materials  subjected  to  shock  pressures  of  less  than  100 
kbar. 

Short  (1965)  proposed  the  reflection  of  stress  waves  from  material  dis- 
continuities as  an  explanation  for  the  occurrence  of  central  peaks.  This 
mechanism  was  based  on  the  partial  reflection  of  the  initial  shock  wave  from 
surfaces  where  the  acoustic  impedance  changes  discontinuously.  These  reflected 
waves,  upon  returning  to  the  crater,  would  reflect  again  from  the  free  surface 
as  tensile  waves,  producing  an  upward  heave  that  would  be  maximum  near  the 
center.  Short  also  suggested  that  if  such  an  effect  did  result  in  central 
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peaks,  then  the  presence  of  such  peaks  on  the  moon  implied  at  least  a zone  of 
higher  acoustical  impedance  below  the  lunar  surface. 

As  an  example,  Short  made  a limited  calculation  of  the  effect  of  this 
mechanism  during  the  impact  that  formed  the  East  Clearwater  Lake  crater  in 
Canada.  He  estimated  that  the  shock  pressure  from  such  an  event  would  atten- 
uate to  a kilobar  at  a depth  of  19  miles  and  if  totally  reflected  would  still 
be  several  hundred  bars  upon  return  to  the  crater  base.  Beals  (1965)  dis- 
counted this  mechanism  because  much  less  than  total  reflection  would  occur. 
Beals  based  his  conclusion  on  the  tables  of  Muskat  and  Meres  (1940),  from 
which  he  inferred  that  for  an  elastic  wave  reflecting  from  the  crust-mantle 
discontinuity  only  about  0.25  percent  of  the  energy  would  be  reflected. 

The  inward  motion  of  material  behind  the  primary  shear  wave  has  also 
been  suggested  as  a cause  of  central  peaks.*  The  motion  behind  that  wave  would 
force  material  into  a smaller  volume  about  the  vertical  axis  through  the 
center  of  the  crater,  because  the  shear  wave  would  be  symmetric  about  that 
axis.  This  squeezing  effect  would  tend  to  cause  upward  velocities  in  a 
manner  similar  to  squeezing  toothpaste  from  a tube.  This  effect  was  evident 
in  calculations  similar  to  the  REVROC  study,  which  showed  upward  motions 
began  where  the  calculated  shear  wave  intersected  the  axis  of  symmetry.  Also, 
as  was  noted  earlier,  Christensen  (1970)  observed  this  mechanism  in  the 
ELK  31  calculation. 

Milton  and  Roddy  (1972)  suggested  that  the  occurrence  of  a central 
peak  in  an  impact  crater  may  indicate  a low-density  impacting  body  such  as  a 
comet.  They  stated  that  a necessary  condition  for  central  peak  formation  may 
be  the  initial  deposition  of  energy  near  the  surface,  and  not  at  some  depth. 
This  condition  may  be  required  because  central  peak  formation  depends  on  a 
complex  interaction  of  the  shock  wave  with  the  free  surface.  If  a major  por- 
tion of  the  initial  energy  is  deposited  too  deeply  into  the  target  material, 
the  region  that  would  have  formed  the  central  peak  would  become  involved  in 
the  crater.  A central  peak  would  form  in  a cometary  impact  crater  because 
the  comet,  consisting  mainly  of  H20  and  C02  ices,  would  volatilize  near  the 
surface  upon  impact,  while  a meteorite  would  penetrate  to  some  depth.  They 
concluded  that  information  on  the  percentage  of  impact  craters  with  central 
peaks  may  indicate  the  ratio  of  large-scale  cometary  impacts  to  meteor  impacts 
on  the  surface  of  a planet. 

*Port,  R.  J. , Air  Force  Weapons  Laboratory,  Personal  Conference,  1973. 
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b.  Gravitational  Mechanisms 

Several  authors  (Shoemaker,  1963;  Dence,  1968;  Gault,  et  al.,  1968) 
have  suggested  that  a cause  of  central  peak  formation  is  a deep  sliding,  or 
base  failure,  resulting  from  the  gravitational  stresses  produced  by  the 
difference  in  height  between  the  rim  and  the  center  of  a crater.  Dence  (1968) 
referred  to  a solution  (figure  7),  which  showed  that  under  a valley  with  walls 


Figure  7.  Theoretical  Lines  of  Slip  under  a Valley  with 
Sides  Dipping  at  30°  to  the  Horizontal  (after 
Dence,  1968) 

sloping  at  30°  the  earth  movement  would  occur  along  slip  lines  that  form  two 
families  of  parabolas  with  the  bottom  of  the  valley  as  the  focus  and  the  apices 
within  the  moving  material.  In  the  cratering  case,  Dence  assumed  the  slip 
lines  would  be  replaced  by  coaxial  surfaces  that  retained  the  upward  turning 
beneath  the  center  of  the  crater.  Motion  along  these  slip  surfaces  would  be 
resisted  by  the  shear  strength  of  the  medium.  There  would,  therefore,  be  a 
minimum  crater  size  for  any  medium  below  which  no  such  motion  could  occur. 

He  described  the  formation  of  a crater  with  a central  mound  as  proceeding 
from  a primary  crater  by  the  walls  sliding  down  and  in  along  deep  slip  sur- 
faces forcing  the  material  under  the  crater  to  bulge. 

Dent  (1974)  has  also  considered  the  failure  mode  of  the  walls  of  a 
crater.  These  modes  of  failure  are  "slope  failure"  (in  which  the  failure 
surface  emerges  in  the  crater  wall)  and  "base  failure"  (in  which  the  failure 
surface  extends  deeply  below  the  bottom  of  the  crater).  He  accomplished  an 
elastic  plane-strain  analysis  of  the  stresses  caused  by  the  excavation  of  a 
semicircular  cavity  at  the  surface  of  a two-dimensional  halfSpace  in  a gravita- 
tional field.  He  concluded  that  with  a Mohr-Coulomb  failure  criterion  the 
slope-failure  mode  would  be  preferred  for  any  size  of  excavation. 
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Pi ke ( 1 971 ) suggested  that  central  peaks  may  be  caused  by  the  centripe- 
tal movement  of  collapsed  rim  material  similar  to  the  Rayleigh  jet  produced  in 
the  transient  craters  in  a liquid  medium.  Harlow  and  Shannon  (1967)  have 
numericaly  simulated  the  splash  of  liquid  drops  into  deep  pools,  showing  the 
development  of  these  splash  jets  in  incompressible  fluids.  Their  results 
showed  that  the  development  of  the  central  jet  was  caused  by  the  gravitational 
collapse  of  the  sides  of  the  crater  into  the  crater  void  provided  that  the 
scale  condition  (gR)1/2/UQ  was  less  than  0.4  where  g is  the  gravitational 
acceleration,  R is  the  radius  of  the  impacting  drop,  and  U0  is  the  impact 
velocity.  However,  these  results  were  changed  significantly  when  compres- 
sibility and  shock  processes  were  involved,  with  the  process  reverting  to 
the  rebo  ind  mechanism  already  described  (Amsden,  1966). 
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Section  III 

NUMERICAL  SIMULATION  OF  SHOCK  WAVE  CRATERING 

The  high-explosive  cratering  experiments  provided  excellent  opportunities 
to  examine  the  causes  of  central  mound  formation  through  numerical  simulation 
because  (1)  the  preshot  material  properties  of  the  medium  were  extensively 
tested,  (2)  the  test  conditions  were  known,  and  (3)  the  post-event  structure  of 
the  craters  and  central  mounds  were  carefully  documented  to  provide  strong 
constraints  on  the  numerical  results.  The  MI, '.ED  COMPANY  II  event  served  as  a 
particularly  useful  experiment  because  of  the  large  size  of  the  central  mound 
compared  to  the  size  of  the  crater.  This  large  size  indicated  that  the  central 
mound  processes  were  particularly  effective  in  this  test  event  and  reduced 
resolution  problems  associated  with  numerical  calculations.  Therefore,  a series 
of  numerical  experiments  was  undertaken  to  simulate  the  MIXED  COMPANY  II  event 
and  determine  the  contribution  of  individual  mechanisms  to  the  formation  of  the 
central  mound.  The  results  showed  that  the  calculation  of  upward  motions  below 
a simulated  crater  was  dependent  on  the  material  compaction  model  in  the  region 
where  strength  effects  were  significant.  The  results  of  one  numerical  simula- 
tion indicated  that  the  presence  of  a lower  "fluid"  material  ma^also  cause 
the  formation  of  a central  mound. 

1.  MIXED  COMPANY  II  EXPERIMENT 

As  previously  stated,  the  MIXED  COMPANY  II  experiment  was  the  detonation 
of  20  tons  of  TNT,  arranged  in  a spherical  charge  of  4 feet  radius,  placed 
above,  and  tangent  to  the  ground  surface.  The  MIXED  COMPANY  test  site  consisted 
of  a thin  deposit  of  sandy  clayey  silt  over  a 70-foot  thick  section  of  Kayenta 
formation  (Ehrgott,  1973).  The  silt,  which  was  1.8  feet  thick  at  the  MIXED 
COMPANY  II  test  site  (Day,  1973)  appeared  to  become  slightly  cemented  at  depth. 
The  Kvyenta  is  a fluvial  deposit  that  consists  of  lenticular  to  irregularly 
bedded  layers  of  fine-to-medium-grained  sandstone,  Siltstone,  and  conglomerate 
with  occasional  layers  or  lenses  of  shale.  The  calculational  models  of  these 
material  were  based  on  properties  determined  from  laboratory  and  CIST  data 
from  the  site*  and  are  given  in  table  3.  This  information  divided  the  test 

*Gajewski , R. , Personal  letter  and  data,  Air  Force  Weapons  Laboratory/DEV, 
Kirtland  AFB,  NM,  1973. 
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Table  3 

MIXED  COMPANY  II  SITE  MODEL* 


Property  Layer 


Depth  to  top  (ft) 

0.0 

1.8 

11.2** 

19.6+ 

Density  (gm/cc) 

1.875 

2.35 

2.47 

2.35 

Compressional  wave  speed  (ft/sec) 

500 

8000 

9000 

8000 

Rarefaction  velocity  (ft/sec) 

1500 

16000 

18000 

16000 

Volume  fraction  of  air  filled  voids 

0.2366 

0.0510 

0.0229 

0.0510 

Poisson's  ratio 

0.25 

0.20 

0.25 

0.20 

Cohesion  (bars) 

0.7 

68 

51 

68 

Angle  of  internal  friction  (°) 

25 

35 

37 

45 

Von  Mises  limit  (kbar) 

0.5 

7.5 

2.1 

11.6 

*Gajewski  R. , Personal  letter  and  data.  Air  Force  Weapons  Laboratory/DEV, 
Kirtland  AFB,  NM,  1973. 

**Changed  to  11.4  ft  in  calculations 

+Changed  to  19.2  ft  in  calculations 

site,  to  a depth  of  60  feet,  into  a 1.8-foot  layer  of  alluvium  over  a halfspace 
that  had  three  layers.  The  properties  indicated  that  the  material  below  the 
soil  had  a much  higher  maximum  yield  strength  than  the  shales  in  the  MIDDLE  GUST 
experiments. 

Observed  crater  morphology  and  structural  information  provided  the  primary 
constraints  on  the  numerical  simulation  of  the  MIXED  COMPANY  II  experiment. 
Detailed  profiles  (figure  8)  of  the  crater  that  was  formed*  showed  that  the 
apparent  crater  extended  a maximum  of  4 feet  below  the  original  ground  level 
at  a radius  from  ground  zero  of  12  feet.  The  crater  was  only  approximately 
symmetric,  with  radii  at  the  original  ground  level  of  20  feet  to  the  north  and 
22  feet  to  the  south.  The  central  mound,  represented  by  true  crater  dimensions, 
extended  a distance  of  8 feet  from  the  vertical  axis  through  ground  zero  and 
was  7 feet  high.  The  mound  was  composed  of  uplifted  and  brecciated  sandstone 
(Roddy,  1973).  A poorly  developed  overturned  flap  and  thin  blanket  of  ejecta 
surrounded  the  crater,  and  no  fused  material  was  found.  Deformation  in  the 


Carnes,  B.  L.,  Personal  letter  and  data.  Waterways  Experiment  Station, 
U S Army  Corps  of  Eng.,  Vicksburg,  MS,  18  October,  1973b. 
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crater  wall  and  rim  consisted  mainly  of  shattering  and  local  brecciation.  A 
piece  of  color-coded  grout,  originally  placed  at  10  feet  depth,  was  found 
1 foot  below  the  top  of  the  central  peak.*  Ground  shock  instrumentation  (Day, 
1973)  included  vertical  and  horizontal  acceleration  gages  at  ranges  of  54  and 
70  feet  and  depths  of  1.5  and  5 feet.  Vertical  and  horizontal  velocity  gages 
were  located  at  similar  depths  and  ranges  greater  than  93  feet.  The  gage  at 
70  feet  range  and  5 feet  depth  indicated  a shock  arrival  at  10  msec.  All  the 
data  indicated  that  a large  signal  was  transmitted  in  the  material  below  the 
depth  of  1.8  feet. 

2.  CALCULATIONAL  MODEL  OF  MIXED  COMPANY  II 

The  numerical  simulation  of  the  MIXED  COMPANY  II  experiment  included  the 
use  of  three  mathematical  models  of  the  physical  processes  that  were  assumed 
to  be  important.  The  first  was  a model  of  the  surface-pressure  boundary  condi- 
tion to  simulate  the  high-explosive  detonation.  The  second  was  a computer  code 
that  modeled  the  initial  response  of  the  ground  to  the  surface  boundary  condi- 
tion. This  code  included  approximations  to  physical  relations  and  the  proper- 
ties of  the  materials  at  the  test  site.  The  final  model  was  a simplified 
ballistic  extension  of  the  conditions  that  were  calculated  using  the  first  two 
models. 

a.  Explosive  Detonation  Model 

The  explosive  detonation  was  modeled  with  data  from  a solution  of  the 
airblast  pressure  as  a function  of  range  and  time  for  a 100-ton  explosive  charge 
on  a rigid  halfspace.**  This  information  was  applied  using  a cube-root  scaling 
procedure  to  provide  a surface-pressure  boundary  condition  for  a 20-ton,  surface- 
tangent  event.  The  procedure  was  to  scale  the  ranges  and  time  of  the  calcula- 
tion by  the  ratio  (100/20) 1/i , apply  the  boundary  condition,  and  tnen  rescale 
range  and  time  by  the  inverse  ratio. 

b.  Ground  Response  Model 

The  initial  response  of  the  ground  was  modeled  with  the  AFT0N-2A  compu- 
ter code  (Niles,  et  al.,  1971).  This  code  models  two-dimensional,  axisymmetric 


♦Personal  conference  with  Major  Lamping  of  the  Air  Force  Weapons  Laboratory 
in  March,  1973. 

**Data  supplied  by  the  Air  Force  Weapons  Laboratory,  see  Port  and  Gajewski, 
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continuum  mechanics  problems  using  elastic-plastic  material  models  that  sim- 
plify to  hydrodynamic  expressions  at  high  internal  energies.  The  theory  of 
this  code  (see  appendix  D)  is  based  on  a specific  method  of  constructing 
finite-difference  approximations  to  the  laws  of  continuum  mechanics  in  inte- 
gral (but  not  necessarily  Lagrangian)  form  that  includes  artificial  viscosity 
(von  Neumann  and  Richtmeyer,  1950)  to  treat  strong  shock  waves.  This  method 
uses  relations  to  describe  mass  conservation,  momentum  conservation,  and  the 
first  law  of  thermodynamics  that  combine  to  describe  also  energy  conservation 
exactly.  AFT0N-2A  is  used  frequently  in  ground-shock  calculations  (see,  for 
example.  Port  and  Gajewski,  1973),  and  the  numerical  errors  associated  with  the 
code  have  been  investigated  extensively  (Cooper,  1971;  Trulio,  et  al . , 1967). 
All  the  calculations  accomplished  during  this  study  used  a Lagrangian  coordi- 
nate system  unless  specifically  stated. 

The  code  provides  information  in  three  forms.  One  form,  termed  a data 
edit,  is  a printed  listing  of  selected  parameters  at  each  calculation  point. 

The  second  form,  called  a restart  dump,  is  a listing  on  magnetic  tape  of  all 
the  information  necessary  to  continue  the  calculation  from  the  time  of  that 
dump.  While  the  primary  purpose  of  this  form  is  to  provide  a restart  capabil- 
ity, these  dumps  also  provide  the  information  required  to  construct  displays, 
termed  flow  field  plots,  of  the  conditions  that  exist  in  the  calculation  space 
at  the  simulated  time  of  each  dump.  The  third  form  is  complete  time  history 
information  of  100  selected  "target"  points.  These  "target"  points  may  be 
considered  to  be  "perfect"  instruments  which  measure  the  forces  and  responses 
of  a mass  particle  without  influencing  the  behavior  of  that  particle.  They 
are  points  that  may  be  located  at  any  position  in  the  calculational  space,  not 
just  at  calculational  meshpoints,  and  move  in  a Lagrangian  manner. 

A description  of  the  calculation  grids  is  required  to  understand  the 
later  ballistic  model  and  information  representations.  The  quantities  in  this 
code  are  computed  on  two  separate,  but  related,  grids.  Motion  quantities  (such 
as  acceleration,  velocity,  and  position)  are  computed  at  the  designated  calcula- 
tion, or  grid,  points.  Thermodynamic  variables  (such  as  stresses,  strains,  and 
internal  energy)  are  computed  at  the  interior  of  the  volume  defined  by  the  four 
surrounding  grid  points.  Thermodynamic  quantities  are,  therefore,  computed  and 
represented  on  a thermodynamic  mesh.  The  combination  of  the  two  grids  divides 
the  volume  surrounding  the  calculation  point  into  four  quarter-volumes  with 
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associated  quarter  masses.  These  quarter-volumes  are  also  used  in  a ballistic 
extension  to  the  code  results. 

The  boundaries  of  the  grid  were  the  surface  boundary,  the  axis  of 
symmetry,  and  two  transmitting  boundaries  (Niles,  et  al.,  1971).  The  trans- 
mitting boundaries  were  imposed  at  60  feet  depth  and  551  feet  range  and  had  no 
significant  influence  on  the  calculated  motions. 

One  of  the  basic  relations  used  to  describe  the  material  properties 
is  the  equation  of  state,  whicn  related  pressure,  P,  to  material  density,  p, 
and  specific  internal  energy.  The  general  equation  of  state  for  the  material 
models  of  the  MIXED  COMPANY  site  was 


P = f(u,  u*) 


(1) 


where  the  excess  compression,  u,  was  defined  as 


P * Ph 


u = 


(2) 


with  p.j  the  initial  material  density  and  u*  the  maximum  excess-compression  ever 
calculated  at  a thermodynamic  mesh  point.  This  functional  relationship  was 
divided  into  a low-pressure  region,  for  u less  than  or  equal  the  volume  fraction 
of  air-filled  voids,  end  a high-pressure  region.  Effects  of  internal  energy 
on  relation  (equation  1)  were  included  by  adding  the  term  Ae  (e  is  the  specific 
internal  energy  and  A is  a constant  assumed  to  be  3 x 10  12  gm/erg)  to  both 
u and  u*.  For  all  calculations,  the  effect  on  pressure  of  this  addition  was 
small.  This  equation  of  state  is  a generalization  of  the  seismic  model 
(Port  and  Gajewski,  1973)  to  allow  for  a permanent  compaction,  as  in  the  CIST 
mode-' . 

The  low-pressure  equation  of  state  was  further  divided  into  a loading 
relation,  for  u equal  u*,  and  an  unloading  relation,  for  u less  than  u*.  The 
loading  relation  was 


P - KL  u 


(3) 
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where 


2 0 1 < 

pi  cp  3 (1  - v! 


defined  for  each  layer  from  the  initial  density,  the  compressional  wave 
speed,  Cp,  and  the  Poisson's  ratio,  v,  of  the  material.  The  unloading  relation 
allowed  for  a linearly  changing  derivative  of  the  equation  of  state  through  the 
relation 


j\  ux  (6  - 1)  + Ku  u + Kv(ux  - u)]  /BUX 


for  (1 


where 


e)  ux  <u  < ux 


u = min  i u I 
X |u,| 


> 3 i/  + i/ 

Ku  Kv 


and  u 3 represented  the  volume  fraction  of  air-filled  voids.  The  parameter 
was  defined  by  the  relation  (equation  4)  where  c was  replaced  by  the  sonic 
velocity  at  the  initial  release  of  pressure,  cu;  and  Ky  was  defined  by  the  same 
relation  with  the  sonic  velocity  as  the  pressure  approaches  zero,  cy.  The 
unloading  hydrostat  was  then 

= Ku[(B  - 1)  ux  + - Kv(u„  - u)»  + KV(BUX)Z 


after  integration  of  equation  5.  For  u less  than  (1  - 3)ux,  the  pressure- 
density  relation  was  assumed  to  be 


p - -kv  [o  - e)  ux  - u] 


with  the  material  in  tension. 

The  low-pressure  region  allowed  for  a reduction  in  the  specific  volume 
of  the  materials  after  a complete  cycle  of  loading  and  unloading  (figure  9). 
The  amount  of  this  reduction  was  defined  by  the  ratio  3 and  the  maximum 


Figure  9.  The  Model  Pressure--Specific.  Volume  Relation  in  the  Low-Pressure  Regime 
for  Layer  2 (Shown  are  the  loading  relation  and  unloading  relations  for 
u*  less  than  U3  and  u*  equal  to  or  greater  than  u3) 
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compression  experienced  within  the  low-pressure  region.  The  parameter 
defined  as 


C = (1  - 3)  (8) 

was  the  compactibil ity  of  a material.  By  transformations  to  specific  volumes, 
the  relation 


(9) 


where  v represents  specific  volume  and  vQ  is  the  zero-pressure  specific  volume 
after  a load-unload  cycle,  allows  the  compactibility  of  a material  model  to  be 
estimated  from  graphs  of  pressure  vs.  specific  volume.  Initial  calculations, 
which  used  the  rarefaction  velocity  from  the  site  data  for  cu  and  cv,  failed 
to  produce  motions  consistent  with  central  mound  formation.  In  the  MIXED 
COMPANY  II  numerical  simulation  these  values  for  the  three  deepest  layers  were 
reduced  to 


cu  = c^  + 1000 

which  implied  a compactibility  of  approximately  10  percent  in  those  layers. 

The  results  of  subsequent  parametric  calculations  showed  the  effects  of  varia- 
tions in  these  two  parameters. 

The  high-pressure  equation  of  state  was  assumed  to  be  independent  of 
u*,  which  resulted  in  one  relation  describing  both  loading  and  unloading.  The 
derivative  of  the  equation  of  state  in  this  region  was  expressed  by 

t ■ - K ' Ku)  ^fnf)  oo) 

where  Km  and  us  were  parameters  determined  from  appropriate  high-pressure  data. 
Relation  10  implied  that 
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P = KL  U3  + yu  - u 3 ) - (Km  - Ku)us 


1 - exp 


/U  3 - U\ 

l us  L 


(11) 


was  the  pressure  density  relationship  for  this  region.  The  values  for  K and 

m 

us  were  determined  from  Hugoniot  data  on  Coconino  sandstone  at  pressures  above 

that  required  to  close  the  air  voids  (figure  10).  The  values  of  K = 680  kbar 

m 

and  u$  = 0.3  provided  the  comparison  shown  in  figure  10  for  layer  2,  with  the 
less  porous  material  having  lower  specific  volumes  for  pressures  below  150  kbar, 
consistent  with  internal  energy  relations.  Although  the  modpl  did  not  compare 
adequately  with  the  data  above  150  kbar,  no  pressure  above  40  kbar  ever 
occurred  in  the  material  during  all  calculations.  The  same  value  of  and 
u$  = u3  + 0.25  were  found  suitable  for  all  four  layers. 

The  shear  modulus,  G,  was  also  calculated  in  tne  equation  of  state 
model.  In  the  low-pressure  region  the  shear  modulus  was  determined  from 


d£ 

du 


'3 

[!  - 2y)  1 

2 

*1+  v) 

(12) 


wnich  allowed  the  Poisson's  ratio  to  remain  constant.  In  the  high-pressure 
region,  the  shear  modulus  was  a constant  defined  by 


G = K 


■(’  - -’'I  i 


ir^j 


(13) 


where  VL  was  the  constant  Poisson's  ratio  of  the  low-pressure  region.  This 
model  is  referred  to  as  the  hybrid  v-G  model  (Zelasko  and  Baladi,  1971; 

Bratton,  1973). 

The  second  basic  relation  used  to  describe  the  materials  is  the  material 
yield  surface.  A simple  Mohr-Coulomb  and  von  Misses  yield  surface  (figure  11) 
that  was  independent  of  the  third  invariant  of  the  deviator  stresses  was  assumed 
for  all  materials.  The  yield  surface,  Y,  was  described  by  the  relation 

tQ  + P tan  4>  P < PYLD 

YMAX  P - PYLD  (14) 
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Figure  11.  Schematic  Yield  Surface  for  Materials 
(Solid  1 ine  indicates  original  yield 
condition  as  a function  of  pressure; 
dashed  line  indicates  the  shifted 
yield  condition;  no  quantitative 
relation  is  expressed) 

with  tQ  the  cohesion,  <p  the  angle  of  internal  friction  the  von  Mises 
yield  strength,  and  PyLD  defined  by 

to  + PYLD  tan  * = YMAX 

Material  separation  was  assumed  to  occur  when  the  value  of  the  yield  surface 
for  a calculation  zone  reached  zero.  At  the  locations  of  material  separation, 
all  forces  except  artificial  viscosity  terms  and  gravity  were  assumed  to  be 
removed.  The  results  of  initial  simulation  attempts  indicated  that,  to  achieve 
motions  consistent  with  the  formation  of  the  observed  crater,  the  yield  descrip- 
tion should  be  shifted  to  remove  the  cohesive  strength  of  the  material  when 
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the  yield  condition  was  first  reached  in  each  thermodynamic  zone.*  This  shift 
was  accomplished  by  the  use  of  a parameter,  S,  evaluated  for  each  zone,  which 
modified  the  expressions  containing  t to 

V1  - S)  (15) 

with  S initially  zero.  The  value  of  S was  incremented  by  0.04  when  the  yield 
condition  was  reached  and  during  the  subsequent  twenty-four  calculational  cycles. 
The  shift  was  accomplished  in  increments  to  avoid  a drastic  change  in  the  yield 
surface  description  during  one  calculational  cycle,  which  might  result  in 
calculational  instability,  and  was  always  completed  in  less  than  0.8  msec  of 
simulated  time. 

Finally,  a relation,  called  a flow  rule,  is  required  to  describe  the 

inelastic  strain  that  occurs  during  flow  with  stress  conditions  limited  by  the 

yield  surface.  The  associated  flow  rule  (Niles,  et  al.,  1971)  was  used  in  the 

calculations  except  as  will  be  noted.  This  flow  rule  was  derived  with  the 

Method  of  Plastic  Potential  (Trulio,  et  al.,  1969)  and  results  in  a plastic 

volumetric  increase,  called  "bulking,"  when  the  yield  surface  is  a function  of 

pressure.  When  the  yield  surface  is  independent  of  pressure,  this  flow  rule 

reduces  to  the  Prandtl -Reuss  flow  rule.  Also,  the  Prandtl-Reuss  flow  rule  was 

used  if  (1)  the  material  was  in  tension,  (2)  the  plastic  volumetric  strain  had 

reached  0.1,  or  (3)  the  value  of  the  yield  surface  was  less  than  0.5  t . The 

o 

first  of  the  conditions  was  caused  by  uncertainties  in  the  description  of  soil 
response  to  tension;  the  second  condition  limited  the  amount  of  bulking;  and 
the  third  condition  was  caused  by  a singularity  in  the  expression  for  the  flow 
rule  when  the  third  invariant  of  the  deviator  stresses  is  ignored  and  the  value 
of  the  yield  surface  is  near  zero. 

The  calculational  grid  spacing  (appendix  A)  was  selected  based  on  the 
decision  that  this  study  was  primarily  interested  in  conditions  in  and  below 
the  crater  region.  Therefore,  the  calculation  grid  and  target  points  were  con- 
centrated in  that  region.  Outside  that  region  the  grid  spacing  was  increased 
geometrically  to  minimize  the  calculation  time.  This  decision  resulted  in  only 
limited  comparisons  between  calculation  and  experiment  instrumentation  data. 


*As  suggested  by  R.  Port,  Air  Force  Weapons  Laboratory,  personal  communication, 
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c.  Ballistic  Extension  Model 

At  the  simulated  time  of  16.4  msec  in  the  MIXED  COMPANY  II  calculation, 
the  material  in  and  below  the  crater  region  was  calculated  to  be  separated  and 
moving  bal 1 i stical ly  (as  will  be  demonstrated).  Because  of  this  complete 
separation,  the  AFTON  calculation  was  stopped  at  16.4  msec  and  a simplified 
ballistic  analysis  was  accomplished  to  estimate  the  final  crater  shape.  This 
analysis  was  accomplished  for  the  region  within  35  feet  range  and  a depth  of 
20  feet  with  the  velocity  conditions  at  16.4  msec  as  initial  conditions.  Each 
grid  point  was  allowed  to  move  ballistically  (appendix  C)  until  the  following 
three  conditions  were  met: 

(1)  The  grid  point  either  immediately  below  or  radially  away  was  not 
moving 

(2)  The  vertical  velocity  was  negative 

(3)  The  density  of  the  material  in  the  bottom,  outward  quarter-volume 
of  the  zone  was  at  least  1.5  gm/cc 

The  motion  of  a grid  point  was  stopped  after  all  these  conditions,  referred  to 
as  the  stopping  criteria,  were  once  met. 

3.  MIXED  COMPANY  II  NUMERICAL  SIMULATION  (MC  2.12) 

The  numerical  simulation  of  the  MIXED  COMPANY  II  experiment,  designated 
MC  2.12,  resulted  in  calculated  flow-field  conditions  at  a simulated  time  of 
8.4  msec  (figures  12  and  13),  which  were  consistent  with  the  formation  of  a 
central  mound.  The  material  within  a range  of  12  feet  and  a depth  of  20  feet 
had  achieved  upward  velocities  with  the  maximum  vertical  velocities  near  the 
vertical  axis.  Also,  all  the  material  within  that  region  had  separated  and 
was  in  ballistic  motion.  A second  velocity  zone,  centered  near  18  feet  range 
and  4 feet  depth,  was  moving  horizontally  outward  and  again  was  completely 
separated.  Only  a flap  of  material  in  the  top  layer  and  beyond  a range  of  18 
feet  had  significant  velocities  and  had  not  separated.  By  16.4  msec  even  this 
flap  was  completely  separated,  with  little  velocity  change  from  the  conditions 
that  existed  at  8.4  msec. 

The  model  crater  (figure  14)  was  formed  by  616.4  msec  and  a fallback  phase 
of  the  problem  was  beginning.  The  radius  and  slope  of  the  model  true  crater 
wall,  defined  by  the  motionless  material  without  extreme  shear  deformation, 
was  consistent  with  the  true  crater  profile.  An  extreme  shear  zone,  with 
horizontal  grid  lines  extending  into  an  overturned  flap,  was  calculated  near  the 
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range  of  the  observed  southern  crater  wall  and  within  the  asymmetry  of  the 
apparent  crater.  The  material  below  the  original  surface  and  between  the  up- 
lift region  and  the  apparent  crater  wall  was  continuing  to  move  ball istical  ly 
toward  the  crater  wall.  Also,  the  bottom  radius  of  the  central  uplift  was 
comparable  to  the  deepest  point  of  the  true  crater,  located  at  a range  of 
8 feet. 

Severa1  deficiencies  of  the  cal culational  results  were  also  apparent. 

First,  the  top  cal culational  line  shown  in  figure  14  at  less  than  30  feet  range 
represented  the  interface  that  was  initially  at  1.8  feet  depth,  while  the  soil 
layer,  even  at  30  feet  range,  was  above  the  2.5  feet  elevation  allowed  in  the 
figure.  However,  the  measurements  of  the  near-surface  motions  in  both  MIDDLE 
GUST  and  MIXED  COMPANY  tests  also  showed  initial  upward  spall  velocities  near 
the  surface  of  at  least  10  ft/sec  caused  by  the  direct-induced  wave  (Bratton, 
1973;  Port  and  Gajewski , 1973).  These  velocities  were  stopped  by  a second  posi- 
tive phase  of  air  overpressure  (Port  and  Gajewski,  1973)  that  would  occur  after 
the  AFTON  calculation  was  stopped  and  was  ignored  in  the  ballistic  extension. 

The  interaction  of  explosion  products,  aerodynamic  forces,  and  particle-particle 
interactions  rendered  a ballistic  treatment  of  ejecta  distribution  irrelevant. 
Also,  the  stopping  of  some  grid  points,  such  as  the  one  located  near  the  18  foot 
range  and  the  0.0-foot  depth,  and  the  continued  motion  of  other  grid  points 
through  the  crater  wall  showed  that  the  stopping  criteria  were  inadequate  for 
model  ejecta. 

Another  discrepancy,  and  of  most  concern  to  a discussion  of  central  peaks, 
was  that  insufficient  upward  motions  appeared  to  have  been  calculated  to  produce 
the  height  of  the  observed  central  mound.  A continuation  of  the  ballistic  cal- 
culation resulted  in  the  highest  grid  point  on  the  symmetry  axis  in  the  figure 
eventually  settling  back  to  2 feet  below  the  original  surface.  As  the  next 
higher  grid  point  was  computed  to  eventually  reach  a height  of  10  feet,  the 
grid  point  in  the  figure  probably  represented  the  top  of  the  calculated  central 
mound  for  this  model.  An  additional  indication  of  the  insufficient  upward 
motion  was  the  small  upward  displacement  of  the  material  at  the  10-foot  depth 
when  compared  to  the  11  foot  upward  displacement  of  the  colored  grout  that  was 
placed  at  that  depth.  Possible  causes  for  this  discrepancy  will  be  discussed 
after  the  results  of  parametric  numerical  experiments  are  described. 

Even  considering  the  displacement  discrepancy,  the  results  were  in  suffi- 
cient agreement  with  observations  of  the  experimental  event  to  warrant  a 
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description  of  the  conditions  that  produced  the  velocity  field  calculated  at 
8.4  msec.  These  conditions  were  shown  by  calculation-space  plots  (figures  15 
and  16)  at  a time  of  3.0  msec.  The  velocity-vector  plot  showed  that  a strong 
clockwise-rotational  velocity  wave  had  reached  a depth  of  14  feet  near  the 
vertical  axis.  Along  the  vertical  axis,  the  motion  changed  from  downward  to 
upward  at  12.6  feet  depth,  or  immediately  behind  the  intersection  of  the 
rotational  motion  with  the  vertical  axis.  Motion  behind  the  wave  was  upward, 
with  inward  and  outward  oscillations  occurring  closer  to  the  origin.  The 
acceleration-vector  plot  showed  clearly  the  wave  structure  in  trie  flow  field. 
The  accelerations  beyond  a radius  of  24  feet  were  directed  radially  outward, 
except  in  the  near-surface  airblast  region,  and  indicated  the  location  of  the 
compressional  wave.  The  acceleration  reversal,  at  a radius  of  24  feet  from 
the  origin,  indicated  the  location  of  maximum  compression.  All  material  within 
that  radius  was  recovering  from  maximum  compressions  and  was,  thus,  described 
by  the  ur^oading  models.  The  wave  that  produced  the  decrease  in  accelerations 
at  the  20-foot  radius  was  caused  by  a second  increase  in  the  overpressure  model. 
The  front  of  the  rotational  velocity  wave  at  16  feet  radius  was  associated  with 
accelerations  which  were  parallel  to  the  wave,  indicating  the  wave  was  the 
principal  shear  wave.  The  velocity  reversal  on  the  vertical  axis  was,  there- 
fore, associated  with  the  location  of  the  principal  shear  wave.  A series  of 
calculated  shear  waves,  produced  by  the  "reflection"  of  the  primary  shear 
wave  at  the  vertical  axis,  extended  toward  ground  zero  and  controlled  the 
sense  of  the  horizontal  velocities. 


These  wave  relationships  were  also  seen  in  the  time  histories  (figures  17 
and  18)  of  the  target  point  located  on  the  vertical  axis,  where  the  geometric 
relations  result  in  the  simplest  analysis  of  vertical  motions.  For  the  target 
point  on  the  vertical  axis  and  at  an  initial  depth  of  10  feet,  the  maximum 
stress  occurred  at  1.4  msec.  The  compressive  wave  was  then  followed  by  a pres- 
sure decrease  that  was  interrupted  at  1.75  msec  by  a combination  of  the  second 
compressive  wave  and  the  reflection  of  the  first  wave  from  the  material  inter- 
face at  11.4  feet  depth.  The  principal  shear  wave  arrived  at  2.3  msec,  pro- 
ducing (1)  a change  in  the  maximum  stress  direction  from  vertical  to  horizontal, 
(2)  a momentary  period  of  a completely  elastic  stress  state,  and  (3)  a reversal 
of  velocity  from  downward  to  upward.  Behind  the  principal  shear  wave,  the 
vertical  stress  was  small  while  horizontal  stresses  were  more  gradually  reduced 
until  the  material  separated  at  5.5  msec.  At  the  time  of  separation,  the 
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Figure  16.  Acceleration  Vector  Plot  for  MC  2.12  at  3.0  msec  after  Detonation 
(Acceleration  is  proportional  to  the  square  of  the  vector  length 
with  the  distance  between  axis  marks  equal  100  kg;  axes  are  in  feet) 
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Motion  Time  History  for  the  Target  Point 
Originally  Located  on  the  Vertical  Axis  at 
10  Ft  Depth  in  the  MC  2.12  Calculation 
(Positive  Y values  indicate  downward  motion) 
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Figure  18.  Stress  Time  History  for  the  Target  Point  of 
Figure  17  with  Compressive  Stresses  Positive 
(Negative  stresses  indicate  material 
separation) 
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target  point  had  an  upward  velocity  of  12  ft/sec,  which  would  result  in  a 
maximum  ballistic  displacement  of  2.2  feet. 

The  horizontal  motions  were  examined  in  the  time  histories  (figures  19  and 
20)  of  the  target  point  initially  located  at  a depth  of  8 feet  and  a distance  of 
4 feet  from  the  vertical  axis.  This  point  was  first  driven  radially  away  from 
the  origin  by  the  compressional  wave.  The  principal  shear  wave,  which  arrived 
near  2.2  msec,  produced  the  inward  velocity  with  a small  enhancement  of  the 
downward  velocity.  The  symmetry  condition  "reflected"  the  shear  wave  at  the 
vertical  axis,  which  resulted  in  a second  shear  arrival  near  2.7  msec  that 
produced  upward  and  outward  motion.  Subsequent  shear  waves  continued  until  the 
material  separated,  at  4.5  msec,  after  which  the  target  point  continued 
ball istically. 

As  noted  before,  a second  region  of  primarily  horizontal  motion  occurred  in 
the  calculation,  centered  near  18  feet  range  and  4 feet  depth.  Time  histories 
in  this  region  (figures  21  and  22)  showed  an  active  period  of  4 msec  after  which 
a condition  of  constant  velocity  was  achieved.  The  constant-velocity  condition 
was  associated  with  no  material  density  decrease,  as  shown  by  the  stresses 
remaining  near  zero  in  contrast  to  the  large  negative  stresses  shown  in  the 
time  histories  below  ground  zero.  Thus,  the  materials  were  moving  radially 
outward,  and  into  a larger  volume,  in  a velocity-range  relation  that  resulted 
in  a constant  density.  This  region  eventually  formed  the  model  true  crater 
wal  1 . 

Two  additional  results  of  this  simulation  were  of  interest  for  shock-wave 
cratering  studies.  First  was  a contour  map  (figure  23)  of  peak  pressures  as 
a function  of  original  position.  This  contour  map  indicated,  by  the  absence 
of  the  40-kbar  contour,  that  the  maximum  pressure  experienced  anywhere  in  the 
calculation!  medium  was  less  than  40  kbar.  The  transmission  of  high  pressures 
in  the  mcd'um  below  the  1.8  foot  soil  layer,  and  the  absence  of  high  pressures 
in  the  soil  layer,  indicated  that  only  direct-induced  signals  were  important 
in  the  lower  medium.  The  increase  of  acoustic  impedance  in  the  layer  between 
11.4  and  19.2  feet  depth  also  resulted  in  modification  of  the  relation  between 
maximum  pressure  and  range.  Second,  the  motion  of  the  target  point  initially 
located  at  a range  of  6 feet  and  a depth  of  1 foot  (figure  24)  simulated  the 
path  of  material  ejected  to  a long  range.  The  particle  path  was  similar  to  the 
motion  of  particles  suggested  by  Gault,  et  al . , (1968)  during  the  excavation 
stage  of  a crater  formed  by  hypervelocity  impact. 
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Figure  19.  Motion  Time  History  of  the  Target  Point  Originally  Located 
at  4 Ft  Range  and  8 .Ft .Depth  in  the  MC  2.12  Calculation 
(Positive  Y values  indicate  downward  motion;  positive  X 
values  indicate  outward  motion;  target  point  path  shown 
with  downward  motion  down  and  outward  motion  to  the 


STRESS  (KBRR) 


al 

0 - 

A - 

itial 

o - 

Stress  Time  History  fo 
Figure  19  with  Compres 
(Negative  stresses  ind 

VEIOL 


AFWL-TR-75-88 


4.  NUMERICAL  PARAMETRIC  STUDY 

The  results  of  the  MIXED  COMPANY  II  numerical  simulation,  MC  2.12,  were 
applicable  only  to  that  experiment.  Extension  of  those  results  to  other 
occurrences  of  central  mounds  required  information  concerning  the  effects  of 
material  properties  on  central  mound  mechanics.  This  information  was  obtained 
through  a parametric  study  of  the  influences  of  compacti bi 1 ity,  layering, 
and  material  yielding  models.  In  addition,  as  a first  attempt  to  apply  the 
study  of  central  mounds  to  an  examination  of  lunar  evolution,  the  influence  of 
a "fluid"  material  below  a solid  layer  was  examined.  The  results  of  the  cal- 
culations in  this  parametric  study  indicated  that  (1)  the  calculation  of  upward 
motions  below  the  crater  was  dependent  on  the  material  compaction  model, 

(2)  the  layering  included  in  the  MIXED  COMPANY  II  simulation  only  slightly 
influenced  the  upward  velocities  below  the  crater,  (3)  the  bulking  model 
included  in  the  associated  flow  rule  contributed  significantly  to  the  calculated 
upward  motions,  (4)  upward  velocities  for  grid  points  on  the  axis  of  symmetry 
were  first  calculated  where  strength  effects  were  important,  and  (5)  the  pres- 
ence of  a lower,  "fluid"  layer  modified  the  calculated  response  in  an  overlying, 
solid  layer  in  a manner  that  ma^  have  eventually  resulted  in  upward  motions. 

a.  Compaction  Model  Effects 

The  results  of  the  DISTANT  PLAIN  6 and  Nuclear  Explosion  numerical 
simulations  had  indicated  that  increased  compaction  of  materials  reduced  the 
calculated  upward  velocities.  The  effects  of  changes  in  the  compaction  model 
on  the  MC  2.12  results  were  examined  in  two  numerical  experiments.  In  the 
first  experiment,  the  values  of  the  initial  unloading  sonic  velocities,  c^, 
in  the  lower  three  layers  were  changed  to  cu  = cl  + 3000  which  increased  the 
compacti bil ity  of  the  materials  to  30  percent.  In  the  second  experiment,  no 
compaction  of  the  material  in  the  lower  three  layers  was  allowed. 

Calculation  MC  2.13,  which  included  the  increased  compacti bil  ity, 
resulted  in  a complete  lack  of  upward  motions  below  the  crater.  By  9.5  msec, 
all  vertical  velocity  components,  (figure  25),  below  the  crater  region  were 
either  downward  or  had  completely  stopped.  Further,  there  were  no  compressive 
stresses  in  the  material  below  the  crater  to  a depth  of  35  feet  to  produce 
significant  upward  velocities  after  this  time.  A ballistic  extrapolation  of 
this  velocity  field  would  be  expected  to  produce  a crater  of  about  22  feet 
radius  and  a maximum  depth  below  the  original  ground  surface  of  3 feet  directly 
below  ground  zero  with  no  central  mound.  These  results  occurred  even  though 
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the  calculated  primary  shear  wave  produced  inward  velocities,  out  not  a change 
in  maximum  stress  direction,  as  it  propagated  down  the  vertical  axis. 

The  results  of  the  calculation  with  incompactibl e lower  layers,  MCP-03, 
indicated  motions  qualitatively  similar  (figure  26)  to  the  MC  2.12  calculation, 
but  with  larger  velocity  magnitudes  below  the  crater.  The  upward  velocity 
after  material  separation  for  the  target  point  on  the  vertical  axis  and 
initially  3 feet  below  the  surface  was  83  ft/sec,  which  was  a factor  of  2 
greater  than  the  velocity  for  the  similar  target  point  in  the  MC  2.12  calcula- 
tion. However,  because  the  magnitude  of  the  upward  velocity  decreased  more 
rapidly  with  depth  in  the  MCP-03  results,  the  velocity  of  the  target  point  at 
the  10-foot  depth  was  only  16.3  ft/sec,  or  a factor  of  1.35  greater  than  the 
equivalent  value  in  the  MC  2.12  calculation.  That  velocity  would  produce  a 
ballistic  displacement  of  only  4.1  feet  upward.  Also,  because  velocity  condi- 
tions in  the  entire  flow  field  were  increased,  the  model  crater  produced  by  a 
ballistic  extension  of  this  calculation  would  be  expected  to  be  less  consistent 
with  the  observed  crater  than  the  MC  2.12  crater  model.  Therefore,  a simple 
reduction  in  the  model  compactibil ity  would  not  provide  an  explanation  of  the 
discrepancy  between  the  heights  of  the  ooserved  and  calculated  central  mound. 

A third  numerical  experiment,  MCP-09.  was  accomplished  because  of  the 
contrast  in  the  results  between  the  MC  2.12  calculation,  with  10  percent  com- 
pactibility,  and  the  MC  2.13  calculation.  The  calculation  MCP-09  was  used  to 
determine  whether  the  upward  velocities  below  the  crater  were  dependent  on  the 
material  compactibil ity  or  the  entire  unloading  relation.  In  the  MCP-09  cal- 


culation, the  same  values  of  cu  as  MC  2.13  were  used,  but  the  compaction  values 


were  reduced  to  the  MC  2.12  values. 


The  results  of  calculation  MCP-09  indicated  that  the  value  of  the  final 
density  after  a cycle  o'-  loading  and  unloading  was  more  important  to  central 
peak  formation  than  an  accurate  description  of  the  unloading  path.  The  calcu- 
lated flow  field  at  9.5  msec  (figure  27)  was  similar  to  the  MC  2.12  calculation, 
with  only  small  differences  existing.  The  differences  included  a maximum 
2v  f>errert  differed  f»  fh**  wirfiffll  Vt'l&ctEte*  jlCng  the  jxK  if 

symmetry.  In  the  MCP-09  calculation,  the  final  velocities  were  all  lower  than 
the  MC  2.12  results  above  9 feet  depth  and  higher  below  that  depth,  which  would 
increase  the  height  of  the  calculated  central  mound.  Also,  the  horizontal 
velocities  were  generally  lower  and  tended  to  be  more  inward  in  the  MCP-09 
calculation. 
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The  similarity  existed  between  the  MC  2.12  and  MCP-09  calculations 
despite  a completely  different  early  history  of  the  two  models.  The  flow 
field  at  2.8  msec  for  MCP-09  (figure  28)  was  much  different  from  the  MC  2.12 
flow  field  (figure  15)  at  3.0  msec.  The  primary  shear  wave  produced  lower  in- 
ward velocities  close  to  the  axis  as  it  passed.  The  maximum  inward  velocity  in 
the  MCP-09  calculation  was  13  ft/sec  at  11.6  ft  depth;  while  in  the  MC  2.12 
calculation,  it  was  28  ft/sec  at  12.3  ft  depth.  Also,  upward  velocities  along 
the  vertical  axis  did  not  occur  below  9 feet  depth,  which  was  significantly 
shallower  than  the  primary  shear  wave  or  even  its  calculated  reflection. 

Finally,  in  contrast  to  the  MC  2.12  calculation,  the  maximum  stress  remained 
the  vertical  stress  after  the  shear  wave  passed.  All  of  these  differences 
were  caused  by  the  different  wave  speeds  of  the  two  models. 

b.  Layering  Effects 

Three  numerical  calculations  were  used  to  determine  the  influence  of 
the  material  layers  that  were  modeled  in  the  MIXED  COMPANY  II  simulation.  In- 
compactible  models  for  the  lower  three  layers  were  used  to  simplify  the  analysis 
of  the  results  and  to  emphasize  che  rebound  effects.  The  first  calculation, 
MCP-01 , consisted  of  a homogeneous,  incompactible  halfspace  with  the  properties 
of  the  second  layer.  In  the  second  calculation,  MCP-02,  the  soil  layer  of 
MC  2.12  was  used  over  the  MCP-01  halfspace.  In  the  third  experiment,  MCP-03, 
the  incompactible  third-layer  model  was  included  from  11.4  feet  to  19.2  feet 
depth,  as  previously  described. 

The  influence  of  the  soil  layer  was  indicated  by  a comparison  of 
velocity-vector  plots  (figures  29  and  30)  near  9 msec.  These  plots  showed 
that  near  the  vertical  axis  the  flow  field  was  similar.  The  flow  was  mostly 
upward,  with  velocities  along  the  vertical  axis  20  percent  greater  for  MCP-02 
than  for  MCP-01.  However,  there  was  much  greater  upward  motion  of  the  region 
centered  near  18  feet  range  and  4 feet  depth  in  MCP-01  than  in  MCP-02.  These 
results  were  caused  by  the  hydrostat  for  the  soil  model.  Since  the  maximum 
pressures  in  the  soil  near  18  feet  were  less  than  100  bars  (figure  23),  the 
soil  responded  as  a soft  material  with  lower  density  than  the  material  below  it. 
Therefore,  the  air  overpressures  were  not  effectively  coupled  at  either  the 
airground  interface  or  the  material  interface  at  this  range.  Furthermore,  the 
soil  layer  was  72  percent  compactible,  which  reduced  the  subsequent  upward 
motion  of  this  layer.  When  the  direct-induced  wave  arrived  in  the  lower  mate- 
rial, a larger  amount  of  the  upward  momentum  was  transferred  to  the  soil  to 
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Figure  30.  Velocity  Vector  Plot  for  MCP-02  at  9.5  msec  after  Detonation  (Velocity 
is  proportional  to  the  vector  length  with  the  distance  between  axis 
marks  equal  100  ft/sec;  axes  are  in  feet) 
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produce  upward  motions  in  the  soil.  However,  at  the  20-kbar  pressures  in  the 
material  below  ground  zero,  both  the  soil  and  the  lower  material  had  similar 
densities  and  bulk  moduli,  which  resulted  in  effective  stress-wave  transmission 
in  this  region.  The  result  was  that  the  soil  layer  only  slightly  affected 
the  strong  direct-induced  signal  in  the  lower  material,  but  reduced  signifi- 
cantly the  effects  of  the  air  shock  away  from  ground  zero. 

Comparisons  betwppn  the  Mf.F-C?  and  KCf  -C3  calculations  showed  that  the 
increased  acoustic  impedance  of  the  third  layer  only  slightly  influenced  the 
upward  motions  that  contributed  to  the  calculated  central  uplift.  The  final 
upward  velocities  fur  target  points  or.  the  vertical  axis  and  initially  located 
from  3 to  6 feet  depth  were  18  ft/sec  higher  for  MCP-03  than  for  MCP-02. 

However,  those  target  points  were  located  in  material  that  was  calculated  to 
be  ejected  from  the  crater.  At  deeper  target  points,  the  relationship  between 
velocities  changed  until  the  final  upward  velocity  of  the  target  point  at 
10  feet  depth  in  MCP-02  was  19  ft/sec,  or  almost  3 ft/sec  higher  than  the 
value  computed  in  MCP-03.  The  final  velocities  for  lower  target  points  on  the 
vertical  axis  were  also  nearly  equivalent  in  the  MCP-02  and  MCP-03  calculations. 

In  the  calculations  that  included  the  third  layer,  the  modeled  von 
Mises  yield  strength  controlled  the  yield  condition  during  maximum  stress  condi- 
tions to  a range  of  6 feet  above  16  feet  depth  in  that  layer.  Although  an 
additional  calculation  indicated  that  this  control  was  not  important  to  the 
calculated  motions  because  of  the  limited  region  involved,  the  change  of  yield 
mode  showed  that  the  calculated  response  would  be  modified  by  increased  stress- 
conditions  . 

c.  Yield  Model  Effects 

The  yield  model  in  the  MC  2.12  calculation  included  a model  of  material 
strength  that  contained  an  assumed  method  of  shifting  the  yield  surface  once 
the  yield  condition  was  met.  This  method  had  not  been  included  in  earlier 
numerical  simulations  of  shock-wave  cratering  events.  Therefore,  the  require- 
ment for  the  shift  and  alternate  ways  of  accomplishing  the  shift  were  examined 
in  three  numerical  experiments. 

The  requirement  for  the  shift  of  the  yield  surface  was  tested  by  a 
calculation,  MC  2.15,  to  examine  if  a complete  absence  of  tensile  strength 
would  also  produce  the  observed  crater  without  shifting  the  yield  surface.  The 
same  model  as  MC  2.12  was  used  in  the  MC  2.15  calculation  except  the  parameter 
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S was  never  changed  from  zero  and  the  material  was  assumed  to  separate  whenever 
negative  pressures  occurred.  This  model  resulted  in  a completely  different 
flow  field  that  was  not  consistent  with  the  observed  crater.  The  plot  of 
velocity  vectors  at  8.4  msec  (figure  31)  showed  that  in  the  second  layer 
the  material  had  large  upward  velocities  only  to  a range  of  8 feet  and  a depth 
of  6 feet.  All  the  upward  velocities  below  that  depth  were  less  than  10  ft/sec. 
Since  there  were  no  calculated  stress  conditions  remaining  at  that  time  which 
would  produce  a change  in  this  flow  field,  the  MC  2.15  model  was  considered  to 
be  inadequate.  This  result  indicated  that  a material  model  with  no  cohesion 
was  required  to  simulate  the  formation  of  the  central  mound  in  the  MIXED 
COMPANY  II  event. 

Alternate  methods  of  accomplishing  the  shift  of  the  yield  surface  were 
tested  in  two  numerical  experiments.  The  MCP-02  calculation  (the  two-layer 
model  of  soil  on  an  incompactible  halfspace)  was  used  as  the  standard  model  for 
these  numerical  experiments.  In  one  calculation,  the  parameter  S was  incre- 
mented by  0.04  only  during  each  of  the  first  25  cal cul ational  cycles  that  the 
yield  condition  was  satisfied.  In  the  second  calculation,  the  parameter  S was 
incremented  by  0.04  only  after  the  calculation  of  the  yield  condition  was  com- 
plete, which  affected  the  yield  surface  only  during  subsequent  calculational 
cycles.  The  calculated  conditions  during  both  of  these  numerical  experiments 
showed  little  variation  from  the  MCP-02  values. 

The  inclusion  of  a yield  model  in  the  computer  code  was  a theoretical 
complication  required  because  materials  have  finite  strengths.  To  demonstrate 
the  effects  of  the  yield  model  on  the  calculated  results,  a no-yield  assump- 
tion was  used  in  one  calculation.  The  MCP-02  model  was  used  as  a standard 
for  this  calculation,  designated  MCP-06,  but  the  yield  condition  in  the  half- 
space below  the  soil  was  completely  ignored. 

A comparison  of  velocity-vector  plots  at  3.2  msec  (figures  32  and  33, 
note  change  in  velocity  scale)  showed  that  the  yield  model  reduced  the  velocities 
experienced  during  ohe  unloading  from  maximum  stresses.  In  the  no-yield  model, 
the  maximum  downward  velocity  and  maximum  normal  stress  of  the  compressional 
wave  near  the  vertical  axis  were  nearly  equal  to  the  analogous  values  when  the 
yield  model  was  included.  However,  in  the  no-yield  case  these  values  were 
immediately  reduced  after  maximum  compression,  so  that  the  downward  velocity 
was  only  5 ft/sec  by  20  feet  depth.  The  change  in  velocity  direction  that 
occurred  at  a 21-foot  radius  from  the  origin  was  associated  with  a change  from 
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Figure  32.  Velocity  Vector  Plot  for  MCP-06  at  3.2  msec  after  Detonation  (Velocity  is 

proportional  to  the  vector  length  with  the  distance  between  axis  marks  equal 
200  ft/sec;  axes  are  in  feet) 
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Figure  33.  Velocity  Vector  Plot  for  MCP-02  at  3.2  msec  after  Detonation  (Velocity  is 

?nnP2rV°na1  t0  the  vector  len9th  with  the  distance  between  axis  marks  equal 
100  ft/sec;  axes  are  in  feet) 
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positive  to  negative  pressures,  which  did  not  occur  in  the  MCP-02  calculation. 
Significant  inward  velocities,  as  high  as  127  ft/sec  at  2.6  feet  depth  and 
17  feet  range,  occurred  before  the  arrival  of  the  shear  wave  caused  clockwise 
rotation.  The  shear  wave  produced  downward  velocities  as  high  as  135  ft/ sec 
near  the  surface  with  inward  velocities  as  high  as  150  ft/sec  at  the  12.3-foot 
depth  and  7-foot  range.  These  values  were  five  times  as  large  as  velocities 
produced  by  the  calculated  shear  wave  in  the  MCP-02  model.  Behind  the  shear 
wave,  the  upward  velocity  on  the  vertical  axis  reached  147  ft/sec  at  12.3  feet 
depth  and  remained  near  100  ft/sec  to  1.8  feet  depth. 

The  yield  model  in  the  MIXED  COMPANY  II  numerical  simulation  included 
the  associated  flow  rule,  the  use  of  which  in  numerical  ground-motion  simula- 
tions is  a recent  development.  As  late  as  1971,  the  Prandtl -Reuss  flow  rule  was 
used  to  describe  plastic  stress-strain  relations  in  numerical  ground-motion 
simulations  (Zelasko  and  Baladi,  1971;  Trulio,  et  al . , 1969).  As  noted  before, 
the  only  difference  between  the  two  flow  rules,  if  the  yield  surface  is  inde- 
pendent of  the  third  invariant  of  the  deviator  stresses,  is  the  bulking  included 
in  the  associated  flow  rule  when  the  yield  surface  is  a function  of  pressure. 

The  effect  of  the  associated  flow  rule  was  examined  in  a numerical  experiment, 
MCP-05,  in  which  the  Prandtl -Reuss  flow  rule  was  used  with  the  MCP-02  model. 

The  MCP-05  calculation  resulted  in  a completely  different  response  of 
the  medium  to  the  applied  boundary  conditions.  Calculated  motion  remained 
predominantly  downward  and  away  from  the  axis  of  symmetry  until  10  msec  with 
no  rotational  motion,  which  was  an  outstanding  feature  of  the  MCP-02  calcula- 
tion at  3.2  msec.  At  10  msec,  inward  velocities  as  high  as  2 ft/sec  were 
calculated  to  a depth  of  26  feet  and  range  of  33  feet;  however,  vertical  motion 
was  still  downward  at  velocities  less  than  10  ft/sec  within  a range  of  15  feet. 
By  the  simulated  time  of  11.4  msec,  the  downward  motion  near  the  vertical  axis 
had  stopped  and  upward  velocities  as  high  as  10  ft/sec  were  calculated  near  the 
vertical  axis  to  a depth  of  20  feet.  Also,  at  this  time,  the  motion  below 
9 feet  depth  was  inward  with  velocities  less  than  2 ft/sec,  and  material  separa- 
tion has  occurred  along  the  axis  to  a depth  of  4 feet.  The  material  separation 
front  proceeded  downward  and  outward  at  a speed  of  3000  ft/sec  until,  by  19.9 
msec,  the  flow  pattern  shown  in  figure  34  has  developed  with  stresses  within 
a 23-foot  radius  from  the  origin  calculated  to  be  near  zero.  During  the  entire 
calculation,  upward  velocities  never  exceeded  15  ft/sec  in  the  halfspace  below 
the  model  soil  layer. 
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Figure  34.  Velocity  Vector  Plot  for  MCP-05  at  19.9  msec  after  Detonation  (Velocity 
is  proportional  to  the  vector  length  with  the  distance  between  axis  marks 
equal  50  ft/sec;  axes  are  in  feet) 
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In  the  MIXED  COMPANY  II  numerical  simulation  and  all  subsequently 
described  calculations,  the  yield  condition  for  the  material  below  the  soil 
layer  was  usually  determined  from  the  Mohr-Coulomb  relation  as  a result  of  the 
high  von  Mises  limits.  Therefore,  the  effect  on  the  MCP-02  calculation  of 
towering  the  von  Mises  liilt  was  tested.  For  this  test,  designated  MCP-£l, 
the  von  Mises  limit  in  the  material  below  1.8  feet  depth  was  set  at  0.21  kbar; 
and,  because  calculated  displacements  were  too  large  to  allow  calculation  with 
rjr^Hr^ly  I vT*rigi*r,  ryt iyst«T  , 'he  generalised  COWdfcate  ejyaUHty 
of  the  AFT0N-2A  code  (Trulio,  1966)  was  used  for  grid  points  initially  at  and 
above  10.6  feet  depth.  The  grid  points  in  that  region  had  fixed  radial  coordi- 
nates but  were  allowed  to  move  vertically  to  retain  proper  stratigraphic  rela- 
tions. The  properties  of  the  mass  that  was  transported  across  calculational 
faces  were  determined  from  the  values  in  the  zone  out  of  which  the  material 
was  being  transported,  and  were  "backward"  transport  terms  (Trulio,  1964; 

Cooper,  1971).  The  coordinate  system  remained  Lagrangian  below  10.6  feet  depth. 

The  time  and  location  of  the  first  upward  motion  on  the  vertical  axis 
were  accurately  determined  in  this  calculation.  Near  the  vertical  axis,  upward 
velocities  were  first  calculated  between  5.6  msec  (figure  35)  and  6.9  msec 
(figure  36)  with  downward  motion  near  the  axis  still  calculated  below  18  feet 
and  above  7 feet  depth.  The  velocity  transition  region,  near  16  feet  depth 
and  5 feet  range,  was  not  associated  with  exceptional  stress  or  acceleration 
conditions  and  only  represented  the  null  region  of  a continuous  decrease  in 
outward  velocities.  Continuation  of  the  calculation  to  10.5  msec  (figure  37) 
showed  that  velocities  with  inward  components  as  large  as  2.5  ft/sec  and  upward 
components  as  large  as  20  ft/sec  were  calculated  below  10  feet  depth.  The 
first  upward  motion  on  the  vertical  axis  occurred  near  12  feet  depth  at  a 
simulated  time  of  6.3  msec. 

The  motion  and  stress  histories  (figures  38  and  39)  of  that  target 
point  were,  therefore,  especially  pertinent  to  an  examination  of  the  mechanics 
of  central  mound  formation.  The  first  signal  was  an  elastic  precursor  that 
arrived  *E  1 l ise:  snd  eauied  the  vtret*  eonlttlyns  to  r.:acl  tkf.  ylelJ  lurfate 
The  plastic  wave  then  arrived  at  1.9  msec,  when  the  maximum  compression  and 
downward  velocity  were  reached.  In  this  plastic  wave  the  stress  deviators  were 
controlled  by  the  von  Mises  yield  condition  and,  therefore,  the  associated  flow 
rule  became  equivalent  to  the  Prandtl -Reuss  flow  rule.  Also,  there  was  a signi- 
ficant difference  between  the  maximum  and  minimum  normal -stress  values.  After 
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igure  36.  Velocity  Vector  Plot  for  MCP-21  at  6.9  msec  after  Detonation  (Velocity 

is  proportional  to  the  vector  length  with  the  distance  between  axis  marks 
equal  100  ft/sec;  axes  are  in  feet) 
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Figure  37.  Velocity  Vector  Plot  for  MCP-21  at  10.5  msec  after  Detonation  (Velocity 
is  proportional  to  the  vector  length  with  the  distance  between  axis  marks 
equal  100  ft/sec;  axes  are  in  ft) 
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the  maximum  stresses  occurred,  the  stresses  were  rapidly  reduced,  under  elastic 
conditions,  interrupted  only  by  the  second  compressional  signal.  A lower 
stress  state  was  then  reached  and  maintained  for  a period  of  3 msec,  during 
which  the  yield  surface  again  controlled  the  deviator  stresses.  This  time, 
however,  the  yield  surface  was  controlled  by  the  Mohr-Coulomb  relation,  which 
resulted  in  volumetric  bulking.  Thus,  while  nearly  constant  stress  conditions 
were  maintained,  because  of  the  downward  momentum  of  the  material  above  this 
target  point,  the  downward  velocity  was  continuously  reduced  until  the  motion 
stopped  at  6.3  msec  and  reversed  direction.  After  6.3  msec,  the  vertical  com- 
pressive stress  was  reduced  while  the  horizontal  normal  stresses  remained  sig- 
nificant, which  resulted  in  elastic  deformation  while  the  vertical  stress 
approached  zero.  The  reduction  of  the  vertical  stress  lowered  the  pressure 
and,  therefore,  the  yield  condition  so  that  the  horizontal  stresses  were 
finally  reduced  by  material  yield  until  material  separation  occurred  at  9 msec. 


The  stress  time  history  of  the  target  point  on  the  vertical  axis  and 
initially  at  20  feet  depth  (figure  40)  provided  an  example  of  the  elastic- 
precursor  development  that  is  actually  observed  in  shock  experiments  (Ahrens 
and  Rosenberg,  1968),  and  served  as  a simple  demonstration  that  strength  effects 
were  properly  treated  by  the  AFT0N-2A  code.  The  first  calculated  stress  wave 
arrived  near  a simulated  time  of  2.5  msec,  after  traveling  most  of  the  distance 
from  the  surface  at  a compressional  wave  speed  of  8000  ft/sec.  The  main  plas- 
tic wave,  however,  travelled  with  the  stresses  controlled  by  the  von  Mises 
yield  condition,  and  the  wave  speed  was 


= (pt) 


1/2 


(16) 


which  from  relation  (4),  became 


Cb  Cp 


1 + v 
3(1  - v) 


1/2 


(17) 


The  wave  speed,  cb,  computed  from  equation  (17)  was  5656  ft/sec,  which  was  con- 
sistent with  the  maximum-pressure  arrival  near  3.5  msec.  Further,  from  Ahrens 
and  Rosenberg  (1968),  the  Hugoniot  elastic  limit  is  related  to  the  von  Mises 
yield  condition  by 
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3GP 

v = e 

YMAX  3B  + 4G  (18) 

with  P the  maximum  elastic  normal  stress  and  B the  bulk  modulus.  Relation 
e 

(18)  was  transformed  to 


P 

e 


2Y 

^tMAX 


1 - v 
1 - 2v 


(19) 


by  elastic  relations  between  B,  G,  and  v.  With  the  values  used  in  MCP-21 , 
relation  ( i 9)  implied  a Pg  of  0.56  kbar,  which  was  the  maximum  normal  stress 
at  the  time  the  von  Mises  yield  condition  was  reached. 

d.  Lower  "Fluid"  Layer 

An  initial  attempt  was  also  made  to  apply  the  study  of  central  mounds 
to  an  examination  of  lunar  evolution.  Two  lunar  craters,  Lansberg  (0°N, 

26. 5°W)  and  Reinhold  (3°N,  23°W)  are  very  similar  except  that  Lansberg  has  a 
central  peak  while  Reinhold  has  no  such  feature.  Reinhold  is  27.5  miles  in 
diameter  and  1.7  miles  deep  (Schmitt,  et  al . , 1967),  and  Lansberg  is  25  miles 
in  diameter  and  1.6  miles  deep  (Eggleton,  1965).  These  two  craters  are  within 
125  miles  of  each  other  (both  southwest  of  Copernicus  in  Oceanus  Porcellarum) 
and  are  probably  located  in  similar  material.  The  main  difference  between  the 
two  craters,  other  than  the  presence  of  a central  mound,  is  that  Lansberg  is 
overlain  by  Oceanus  Porcellarum  material  while  Reinhold  overlies  that  material 
(Schmitt,  et  al . , 1967). 

One  explanation  of  the  relationship  between  these  two  craters  is  that 
conditions  favoring  central  mound  formation  existed  at  the  time  Lansberg  was 
formed  and  were  sufficiently  modified  by  the  time  of  the  Reinhold  event  to 
prevent  the  formation  of  a similar  mound.  Perhaps,  at  the  time  Lansberg  was 
formed,  a solid  crust,  approximately  one  crater  radius  thick,  existed  over  a 
molten  layer;  while,  by  the  time  Reinhold  was  formed,  the  crust  was  much  thicker. 
Such  a relation  of  a growing  crust  is  contained  in  several  published  thermal 
models  (Toksoz,  et  al.,  1972;  McConnell  and  Gast,  1972)  and  is  used  by  Simmons, 
et  al.,  1973,  to  account  for  the  discontinuous  increase  in  seismic  velocity  at  a 
depth  of  15  miles  with  a constant  seismic  velocity  in  the  depth  interval  15  to 
30  miles  on  the  moon. 
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The  impact  calculators  required  to  test  fully  such  a model  of  the 
relation  between  these  two  craters  were  not  accomplished.  However,  a calcula- 
tion, MCP-12,  was  completed  that  inserted  a model  “fluid"  halfspace  below 
16.2  feet  depth  in  the  MCP-02  model  to  provide  a preliminary  test  of  the  in- 
fluence of  such  a medium.  In  the  MCP-12  calculation,  the  "fluid"  was  modeled 
by  an  incompatible  hydrostat,  with  the  values 

Initial  density  - 2.47  gm/cc 
Seismic  velocity  - 6000  ft/sec 
u 3 - 0.0229 
Poisson's  ratio  - 0.45 

and  to  simulate  viscosity  a constant  yield  surface  value  of  34.5  bars  for  com- 
pressive pressures  with  material  separation  assumed  for  negative  stresses. 

Also,  since  the  transmitting  boundary  assumed  impinging  elastic  waves,  the 
bottom  boundary  was  located  at  156.2  feet  depth.  A special  "sliding"  interface 
condition,  which  limited  the  shear  stress  transmitted  across  the  interface  to 
the  shear  capability  of  the  fluid  (Niles,  et  al.,  1971),  was  used  at  the  solid- 
fluid  interface. 

The  calculated  flow  field  at  16.6  msec  (figure  41)  showed  that  the 
"fluid"  layer  extensively  influenced  the  motions  of  the  solid  material  above 
the  interface.  This  effect  was  particularly  strong  in  a cylindrical  region  of 
15  feet  range  and  below  9 feet  depth  where  the  motion,  instead  of  being  upward 
as  in  the  MCP-02  calculation,  was  downward  with  velocities  as  high  as  70  ft/sec. 
This  motion  was  producing  large  downward  material  displacements  with  the  cylin- 
der of  solid  material,  which  was  completely  separated,  moving  into  a developing 
depression  in  the  "fluid."  The  motion  in  this  "fluid"  was  similar  to  that  cal- 
culated by  Harlow  and  Shannon  (1967)  for  the  splash  of  a liquid  drop  into  a deep 
pool,  and  an  extension  of  this  calculation  would  be  expected  to  develop  large 
upward  velocities  near  the  axis  of  symmetry  as  the  fluid  recovered  under  gravi- 
tational flow.  The  amount  of  this  recovery  would  be  controlled  by  the  properties 
of  the  fluid,  which  is  speculated  to  control  also  the  development  of  a central 
mound  formed  by  the  solid  material  as  that  material  is  redirected  upward  by  the 
fluid. 
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41.  Velocity  Vector  Plot  for  MCP-12  at  6.6  msec  after  Detonation  (Velocity 
is  proportional  to  the  vector  length  with  distance  between  axis  marks 
equal  100  ft/sec;  axes  are  in  feet) 
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SECTION  IV 

DISCUSSION  AND  CONCLUSIONS 


i 


1.  DISCUSSION 

a.  Evaluation  of  Results 

Numerical  modeling  of  physical  events  is  subject  to  two  basic  causes  of 
error.  One  of  these  causes  is  the  numerical  error  caused  by  replacing  the  space 
and  time  continuums  with  a discrete  grid  on  which  calculations  are  accomplished 
at  specific  moments  in  time.  This  form  affects  the  accuracy  of  the  calculated 
numerical  values  and  usually  can  be  reduced  by  decreasing  the  grid  spacing  and 
cal culational  time  increments.  Evaluation  of  this  form  of  error  in  complex 
calculations  can  only  be  inexact;  however,  from  previous  investigations  of 
numerical  error  in  the  AFT0N-2A  code  (Cooper,  1971;  Trulio,  et  al.,  1967),  such 
errors  are  not  expected  to  seriously  affect  relations  between  calculated  values. 
The  second,  and  more  serious,  cause  of  error  is  the  use  of  invalid  mathematical 
models  for  the  actual  physical  processes.  This  cause  of  error  can  result  not 
only  in  numerical  values  that  are  incorrect,  but  also  in  a calculated  response 
that  is  completely  invalid. 

The  MIXED  COMPANY  II  numerical  simulation,  MC  2.12,  was  subject  to 
modeling  errors  of  both  forms  in  the  three  models.  The  first  was  the  model  of 
the  overpressure  boundary  condition  used  to  simulate  the  high-explosive  event. 
This  model  has  been  suggested  (Trulio  and  Perl,  1977)  to  underestimate  the  over- 
pressure impulse  of  the  MIDDLE  GUST  III  event  by  at  least  40  percent  at  ground 
zero.  The  second  was  the  model  of  the  ground  response  to  the  overpressure 
boundary  conditions  during  the  first  16.4  msec.  This  model  included  only 
approximations  of  (1)  the  test  site,  (2)  the  properties  of  the  materials  at 
that  site,  and  (3)  descriptions  of  physical  relations.  The  third  model  was  the 
ballistic  extension  of  the  calculated  conditions  at  16.4  msec. 


A necessary  but  not  sufficient  condition  for  the  models  to  provide  a 
simulation  of  the  actual  physical  processes  is  that  the  calculated  final  condi- 
tions are  the  same  as  the  final  conditions  observed  after  the  physical  event. 
The  results  of  the  MIXED  COMPANY  II  numerical  simulation  do  not  completely  meet 


this  condition;  however,  the  calculation 


may  still  be  useful  to  a study  of 
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central  peak  mechanics  if  the  inadequacies  of  the  models  do  not  include  the 
causes  of  central  mound  formation. 

Thus,  the  interpretation  that  the  model  produced  insufficient  displace- 
ments below  ground  zero  is  considered  to  be  the  major  discrepancy  of  the  simula- 
tion. Since  the  ballistic  extension  model  is  not  expected  to  underestimate  the 
displacements,  the  cause  of  this  discrepancy  is  thought  to  exist  in  at  least 
one  of  the  first  two  models.  One  possible  cause  is  that  the  impulse  of  the 
model  overpressure  was  too  low.  An  increase  of  the  overpressure  impulse  would 
cause  a stronger  shear  wave  and  more  bulking,  because  the  stresses  and  the  com- 
pressional  pressure  duration  would  be  greater.  Another  possible  cause  is  that 
the  dependence  of  the  yield  surface  on  the  third  invariant  of  the  deviator 
stresses  was  not  included  in  the  yield  model.  It.  many  natural  materials,  the 
yield  surface  is  dependent  on  the  third  invariant  (White,  1973),  which  indicates 
that  the  yield  model  used  in  the  numerical  calculations  was  incomplete.  Other 
possible  causes  of  the  discrepancy  include  errors  in  the  values  used  to  model 
the  test  site  and  numerical  error. 

b.  Mechanical  Model  of  Central  Mound  Formation 

However,  the  development  of  a model  of  the  causes  of  central  mound 
formation  is  warranted  because  the  final  conditions  calculated  in  the  MIXED 
COMPANY  II  numerical  simulation  are  sufficiently  similar  to  the  observed  condi- 
tions after  the  physical  experiment.  The  results  of  the  numerical  experiments 
(table  4)  can  be  used  as  a guide  to  this  model.  In  this  model,  the  principal 
cause  of  central  mound  formation  is  the  rebound  of  material  following  the  maxi- 
mum shock  pressure.  If  this  rebound  is  sufficiently  reduced  by  compaction  of 
the  materials,  a central  mound  will  not  form  in  the  crater  unless  there  is  a 
material  below  the  bottom  of  the  crater  that  responds  like  a fluid  to  produce  a 
gravitational  splash  jet.  The  rebound  must  occur  in  the  region  where  material 
strengths  are  sufficient  to  retard  some  of  the  outward  displacements  caused  by 
the  initial  compressional  wave.  The  formation  of  the  central  mound  is  enhanced 
by  the  bulking  of  material  as  the  material  brecciates  and  by  the  inward  motion 
caused  by  the  principal  shear  wave.  Inward  displacements,  however,  are  associ- 
ated with  a ductile  response  of  the  material.  The  presence  of  material  dis- 
continuities below  the  forming  crater  will  reflect  stress  waves  that  might  add 
to  the  formation  of  a central  mound.  Because  of  tne  reduced  density  in  the 
region  forming  the  central  uplift,  a period  of  material  separation  exists  that 
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Table  4 

SUMMARY  OF  NUMERICAL  SIMULATION  RESULTS 


Model 


Basel ine 


Designation  Mode] 
MC  2.12*  Numer 


Model 

Change 


Numerical  simulation  of  MIXED 
COMPANY  II  with  a hydrostatic 
compaction  factor  (page  41)  of 
10  percent 


Major 
Resul  t 

Generally 
consistent  with 
observations 


2. 

MC  2.13* 

MC  2.12 

30%  compaction 

No  central  mound 

80 

3. 

MCP-09* 

MC  2.13 
MC  2.12 

10%  compaction 
different  release 
adiabat 

Like  MC  2.12 

84 

4. 

MC  2.15 

MC  2.12 

Cohesion  retained 

No  central  mound 

91 

5. 

MCP-02* 

No  compaction 
two  layers 

- 

Central  mound 

87 

6. 

MCP-01 

MCP-02 

No  soil 

Little  change  in 
central  mound 

87 

7. 

MCP-03 

MC  2.12 
MCP-02 

No  compaction 
4 layers 

Increased 

velocities 

Little  change 

83 

90 

8. 

MCP-05* 

MCP-02 

No  bulking 
(Page  45) 

Completely 

different 

97 

9. 

MCP-06 

MCP-02 

No  yield 

Large  transient 
velocities 

9G 

10. 

: 

MCP-07 

MCP-02 

Shift  during 
yield 

No  change 

93 

.11. 

MCP-08 

MCP-02 

Shift  after 
yield 

No  change 

93 

12. 

MCP-12* 

MCP-02 

"Fluid"  layer 

Probable  fluid 
splash 

112 

13. 

i 

MCP-21 * 

MCP-02 

0. 21  kbar  yield 

Deeper  upward 

99 

important  implications  for  central  peak  formation 


motions 


may  allow  the  gravitational  sliding  of  the  crater  walls  along  deep  slip  surfaces 
to  produce  additional  inward  displacements. 

c.  Applications  to  Previous  Information 

The  general  applicability  of  this  mechanical  model  of  central  mound 
formation  can  be  tested  by  comparisons  with  previous  observations  of  central 
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mound  occurrence.  Such  a model  should  also  provide  a means  of  understanding 
and  evaluating  previous  calculations. 


Comparisons  should  first  be  made  with  the  high-explosive  experiments 
to  determine  if  such  a model  is  successful  under  loading  conditions  similar  to 
those  for  which  the  model  was  first  determined.  The  absence  of  central  peaks 
in  the  20-ton  detonations  in  Canada  would  be  attributed  to  the  compactibil ity 
of  the  top  25  feet  of  soil  at  that  site.  The  presence  of  central  mounds  in  the 
higher-yield  experiments  would  indicate  that  the  lower,  less-compactible  mate- 
rial was  being  significantly  influenced.  The  ductile  characteristics  of  the 
deformation  in  those  central  mounds  can  either  be  attributed  to  the  independence 
of  the  yield  surface  on  confining  pressure  or  to  a "fluid"  action  of  the  mate- 
rial below  25  feet  depth.  The  occurrence  of  a central  mound  in  the  100-ton 
spherical  event  and  not  in  the  100-ton  hemispherical  event  can  be  attributed 
to  an  increased  precompaction  of  the  soil  below  the  spherical  charge  combined 
with  an  increased  overpressure  impulse  at  ground  zero.  The  mechanical  model  is 
also  applicable  to  the  MIXED  COMPANY  test  series.  The  MIXED  COMPANY  I and  III 
high-explosive  detonations  certainly  caused  higher  stresses  in  the  medium 
(which  caused  more  ductile-like  behavior  of  the  layer  at  12  feet)  than  the 
MIXED  COMPANY  II  experiment.  This  ductile-like  behavior  resulted  in  less- 
pronounced  central  mounds  in  the  first  and  third  experiments  than  the  large 
breccia  cone  of  the  MIXED  COMPANY  II  event.  Further  evidence  of  the  effect  of 
a change  in  yield  mode  is  provided  by  the  central  mound  of  the  MIXED  COMPANY  I 
experiment  with  the  intact  flanks  and  the  brecciated  central  core.  The  MIDDLE 
GUST  series  is  less  easily  interpreted  because  of  the  jointing  in  the  test 
sites;  however,  the  test  series  does  indicate  that  the  more  ductile-like  behav- 
ior of  the  material  at  the  MIDDLE  GUST  sites  than  at  the  MIXED  COMPANY  site 
resulted  in  more  subdued  central  mounds.  The  transition  from  central  mounds 
to  the  central  trough  of  the  MIDDLE  GUST  III  event  was  caused  by  the  higher 
stresses  that  reached  the  lower  shale  layer.  Comparisons  between  the  three 
test  series  illustrate  the  influence  of  compactibil  ity  on  the  occurrence  of 
central  mounds  and  of  yield  mode  on  the  deformation  in  the  central  mound. 

Comparisons  with  impact  structures  provide  a test  of  the  applicability 
of  the  model  to  a different  form  of  loading.  As  the  presence  of  shatter  cones 
indicates,  the  maximum  stresses  in  the  central  mounds  were  above  but  near  the 
Hugoniot  elastic  limit,  producing  ductile  failure  with  eventual  inward  displace- 
ments. The  change  in  the  yield  mode  at  increased  depths  results  in  a central 
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mound  with  a breccia  core  of  reduced  density.  The  dependence  on  gravity, 
suggested  by  Hartmann  (1972),  would  be  caused  by  the  initial  overburden 
stresses  reducing  the  porosity  of  deeper  material,  combined  with  the  additional 
rebound  caused  by  the  dynamic  excavation  of  the  crater.  The  increased  pressures 
caused  by  the  lithostatic  load  would  also  cause  a higher  Mohr-Coulomb  yield 
value,  which  would  increase  the  bulking  and  shear-wave  magnitude.  The  change 
from  central  peaks  to  peak  rings  may  be  caused  by  a gravitational  collapse  of 
the  breccia  core  or  a complete  penetration  of  a solid  crust.  The  transition 
from  simple  to  complex  impact  structures  in  Canada  may  be  caused  by  similar 
gravity  influences  combined  with  a weathering  of  the  Canadian  shield  to  a 
depth  of  one  to  two  miles,  with  an  increased  jointing,  porosity,  and  a reduced 
material  strength. 

A further  test  of  the  model  is  actually  provided  by  its  use  to  explain 
why  central  peaks  do  not  occur  in  some  shock-wave  cratering  events.  The 
Barringer  crater  has  no  central  mound  because  it  was  formed  in  porous  sandstone. 
All  the  nuclear  tests  occurred  in  or  over  compactible  material.  The  hyper- 
velocity impact  experiments  in  loose  sand  targets  produced  only  simple  crater 
shapes . 

An  understanding  of  the  results  of  previous  calculations  is  facilitated 
by  use  of  the  mechanical  model.  Since  the  ELK-31  soil  model  was  precompacted 
before  the  beginning  of  the  calculation,  upward  motions  were  calculated  only 
in  that  numerical  simulation  of  the  DISTANT  PLAIN  6 experiment.  Upward  motions 
were  calculated  in  the  Sierra  Madera  numerical  simulation,  which  was  considered 
to  successfully  model  the  event.  However,  the  motions  in  this  calculation 
resulted  from  invalid  physical  models.  In  particular,  a no-strength  condition 
was  assumed  once  the  calculated  pressure  in  a zone  was  negative.  Because  of 
this  assumption,  the  calculation  actually  simulated  the  response  to  gravita- 
tional force  of  a perfect  fluid  surrounding  an  initial  void.  Evidence  sup- 
porting this  interpretation  of  the  calculation  is  that  the  calculated  motions 
required  a simulated  time  of  15  seconds  to  develop  fully  and  were  still 
prevalent  at  a simulated  time  of  30  seconds.  Errors  in  the  calcul ational  model 
which  probably  reduced  the  earlier  upward  velocities  included  the  use  of  in- 
formation from  only  a Hugoniot  curve  to  determine  all  the  parameters  of  an 
equation  of  state,  and  an  initial  yield  strength  of  0.2  kbar  that  was  indepen- 
dent of  pressure. 
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2.  CONCLUSIONS 

Based  on  results  of  the  numerical  calculations,  four  major  conclusions  were 
reached: 

a.  A central  peak  will  form  in  a crater  in  a solid  medium  when  the  rebound 
of  material  from  maximum  compression  is  sufficient. 

b.  The  rebound  mound  must  occur  in  the  region  of  the  medium  where  material 
yields  but  the  material  strength  is  significant  compared  to  the  maximum 

stresses. 

c.  The  rebound  mound  is  significantly  enhanced  by  plastic  volumetric 
increases  that  occur  during  yielding  when  material  strength  is  a func- 
tion of  pressure. 

d.  The  presence  of  a fluid-like  material  below  a solid  layer  nnght  also 
cause  the  formation  of  a central  peak  in  a shock-produced  crater. 

These  conclusions  form  the  basis  of  a model  of  the  mechanics  of  central  peak 
formation. 

Because  the  observed  occurrence  and  structural  relations  of  central  peaks 
in  craters  produced  by  high-explosive,  nuclear,  and  hypervelocity  impact  events 
can  be  explained  using  the  mechanical  model,  I conclude  that  this  model  de- 
scribes the  causes  of  central  mound  formation  in  shock-wave  cratering  events. 
These  peaks  are  structures  that  are  controlled  primarily  by  the  properties  of 
the  cratered  medium.  The  properties  of  the  source  of  energy  for  the  cratering 
event  influence  central  mound  formation  only  by  modifying  the  stress  distribu- 
tion in  the  cratered  medium.  Shear  waves  and  stress-wave  reflections  from 
discontinuous  increases  in  acoustic  impedance  enhance  central  mound  formation 
but  are  not  the  primary  cause.  Deep  gravitational  sliding  may  occur  and  in- 
crease the  inward  displacement  of  material,  but  this  process  is  an  effect 
rather  than  a cause  of  central  mound  formation. 

Also,  since  the  occurrence  and  structure  of  a central  uplift  are  controlled 
by  the  properties  of  the  cratered  medium  at  the  time  of  the  shock-wave  cratering 
event,  these  structures  provide  a record  of  those  properties.  With  the  mechani- 
cal model  as  a guide,  this  record  may  be  useful  in  examining  the  evolutionary 
and  spatial  relations  of  the  media  in  which  shock-wave  cratering  events  have 
occurred.  Evolutionary  models  of  solar-system  bodies  can  be  constrained  by  the 
structure  of  craters  observed  in  spacecraft  photography. 
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Because  the  processes  of  the  mechanical  model  are  described  in  the  calcu- 
lational  code,  that  code  becomes  a tool  that  can  be  used  to  examine  occurrences 
of  central  peaks.  Such  an  application  to  larger-scale  events  with  higher  maxi- 
mum pressures  and  different  loading  conditions  will  require  extensive  modifica- 
tions to  the  actual  code  used  during  this  study.  These  changes  include  a more 
general  equation  of  state,  an  active  inclusion  of  initial  lithostatic  stresses, 
and  a different  scheme  of  grid  definition.  However,  the  modifications  can  be 
made  within  the  present  calculational  structure  of  the  AFT0N-2A  code. 

3.  FUTURE  WORK 

The  results  of  the  MCP-12  numerical  experiment  indicate  that  the  presence 
of  a "fluid"  below  a solid  crust  will  significantly  change  the  response  of  the 
crust  to  a shock-wave  cratering  event.  The  manner  of  this  change  is  speculated 
to  eventually  result  in  the  development  of  upward  motions  in  the  solid  medium. 
Therefore,  the  structural  relations  between  the  craters  Lansberg  and  Reinhold 
might  be  the  result  of  a growing  lunar  crust.  A program  to  simulate  numerically 
the  events  which  caused  these  two  craters  is  recommended.  While  such  a simula- 
tion could  be  completely  accomplished  with  the  AFT0N-2A  code,  it  is  recommended 
that  only  the  early  part  of  the  Lansberg  calculation  be  accomplished  with  that 
code  to  develop  the  velocity  distribution  in  the  lower  layer.  The  results  of 
this  part  would  then  be  used  as  initial  conditions  in  a code  that  explicitly  and 
more  efficiently  treated  viscous  flow. 

Both  laboratory  and  numerical  experiments  should  be  accomplished  to  deter- 
mine if  the  momentum  and  energy  density  of  a hypervelocity  impact  have  any 
special  influence  on  the  general  structural  characteristics  of  an  impact  crater. 
Since  crater  size  has  been  related  to  the  energy  of  the  event  (Baldwin,  1973), 
if  the  influence  of  momentum  could  be  determined  then  the  mass  and  velocity  of 
the  impacting  body  could  be  estimated.  If  the  influence  of  energy  density  on 
crater  structure  could  be  determined,  then  the  density  of  the  impacting  body 
could  also  be  determined. 

In  the  yield  model  used  to  accomplish  the  calculations,  the  assumption  that 
the  yield  surface  was  independent  of  the  third  invariant  of  the  deviator 
stresses  is  not  generally  valid.  The  general  associated-flow-rule  is  dependent 
on  the  third  invariant  (Trulio,  et  al.,  1969),  and  at  least  one  model  of  the 
dependence  of  the  yield  surface  on  the  third  invariant  exists  (White,  1973). 

The  need  for  including  this  addition  in  the  AFT0N-2A  yield  model  should  be 
examined  further. 
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Data  on  occurrence  of  central  peaks  should  be  reviewed  for  spatial  and  age 
relationships.  For  example,  if  the  lunar  craters  of  25  miles  diameter  that 
had  central  mounds  has  similar  age  characteristics,  then  conditions  favorable 
to  central  mound  formation  would  have  occurred  at  the  time  those  craters  were 
formed.  Then,  if  that  time  varied  with  crater  size,  the  evolution  of  those 
conditions  could  be  examined.  Also,  Manley,  et  al.,  1973,  reported  that  the 
distribution  of  central  peaks  in  the  cratered  areas  on  Mars  is  characterized 
by  clusters  and  seems  nonrandom.  They  stated  that  the  curve  of  central  peak 
frequency  versus  crater  diameter  has  a maximum  at  a crater  diameter  of  5 to 
7.5  miles  and  decreases  rapidly  as  the  diameter  increases.  They  also  reported 
that  large  numbers  of  central  peaks  occur  in  craters  of  2.5  to  7.5  mile 
diameters,  with  10  percent  of  those  craters  having  multiple  peaks.  The  physical 
model  of  Mars  should  be  examined  for  special  conditions  which  would  favor  cen- 
tral mound  formation  between  1 and  4 miles  deep.  Similar  examinations  should 
also  be  made  of  the  photographic  information  on  Mercury.  Such  studies  would 
provide  important  guides  to  evolutionary  models. 
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APPENDIX  A 

COMPUTATION  PARAMETERS 


The  complete  set  of  computation  parameters  for  all  the  numerical  experiments 
described  in  the  text  may  be  divided  into  two  categories.  The  first  category 
includes  the  parameters  that  were  not  varied  during  any  of  the  problems.  These 
parameters  are  given  once  in  this  appendix  and  apply  to  all  the  numerical 
experiments.  The  second  category  includes  the  parameters  that  were  varied  for 
at  least  one  of  the  numerical  experiments.  These  parameters  are  given  for  the 
reference  calculations  and  the  variations  are  described  for  each  of  the  other 
numerical  experiments. 

CONSTANT  NUMERICAL  PARAMETERS 

1.  Explosive  yield  - 20  tons  TNT 

2.  Initial  time  - 0.206  msec 

3.  Gravitational  acceleration  - 32.2  ft/sec 

4.  Artificial  viscosity  coefficients 

Bulk  linear  constant  - 0.06 
Bulk  quadratic  constant  - 4.0 
Deviatoric  quadratic  constant  - 4.0 
Deviatoric  linear  constant  - 0.06 

5.  Calculation  grid  definition 


RADIAL  COORDINATES 

Initial 

(ft) 

Number 

Initial 
Spacing 
111) 

Growth 

Rate 

First  Region 

0 

21 

1.0 

1 .0 
1.1 

Second  Region 

20 

40 

1.2 

VERTICAL 

Initial 

m 

COORDINATES 

Initial 

Spacing 

Number  (ft) 

Growth 

Rate 

First  Region 

0. 

4 

-0.6 

1 .0 

Second  Region 

-1.8 

22 

-0.8 

1 . 0 

Third  Region 

-19.4 

16* 

-1.13 

1 . 1 

*For  MCP-12,  changed  to  27. 
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6.  Initial  target  point  coordinates  - 100  locations  formed  by  combination  of 
X Coordinate  (ft)  0.,  2.,  4.,  6.  , 8. , 10.,  14.,  18.,  22.,  26. 


Y Coordinate  (ft)  -1., 

-2. , -3. , -4. 

, -6. , -8. , 

-10.,  -12. 

, -16.,  -20 

Material  properties 

Layer  1 

Layer  2 

Layer  3 

Layer  4 

Initial  density 
(gm/cc) 

1.875 

2.35 

2.47 

2.35 

u3 

0.2366 

0.051 

0.0229 

0.051 

us 

0.50 

0.30 

0.27 

0.30 

Km(kbars) 

680 

680 

680 

680 

Sublimation  energy 
(xlO12  ergs/gm) 

0.02 

0.02 

0.02 

0.02 

A (xl 0~ 1 2 gms/erg) 

3 

3 

3 

3 

8.  Problem  specifications 

MC  2.12 


Layer  1 

Layer  2 

Layer  3 

Layer  4 

Depth  to  top  (ft) 

0.0 

-1 .8 

-11.4 

-19.4 

cL(ft/sec) 

500 

8000 

9000 

8000 

cu (ft/sec ) 

1250 

9000 

10000 

9000 

cv(ft/sec) 

500 

8000 

9000 

8000 

Poisson's  ratio 

0.25 

0.20 

0.25 

0.20 

tQ(bars) 

0.7 

68 

51 

68 

tan  4> 

0.466 

0.7 

0.75 

1.0 

VMAX  <kbars> 

0.5 

7.5 

2.1 

11.6 

Tension  limit 
(bars) 

0 

-68 

-68 

-68 

MC  2.13 

: MC  2.12  with  the  changes 

Layer  1 

Layer  2 

Layer  3 

Layer  4 

cu(ft/sec) 

1250 

11000 

12000 

11000 

2.15:  MC  2.12  with 

no  yield  surface  shift 

and  the 

changes 

Layer  1 

Layer  2 

Layer  3 

Layer  4 

Tension  limit 

0 

0 

0 

0 

(bars) 

100 
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MCP-05 : 
MCP-06: 
MCP-07 : 
MCP-08: 
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MCP-01:  MCP-02  with  Layer  2 in  the  Layer  1 position  also 

MCP-02 


Layer  1 

Layer  2 

Depth  of  top  (ft) 

0.0 

-1.8 

cL(ft/sec) 

500 

8000 

cu( ft/sec) 

1250 

8000 

cy( ft/sec) 

500 

8000 

Poisson's  ratio 

0.25 

0.20 

tQ(bars) 

0.7 

68 

tan  <p 

0.466 

0.7 

ymax  (^ars ) 

0.5 

7.5 

Tension  limit  (bars) 

0 

-68 

MCP-03:  MCP-02  with  the  changes 


Depth  to  top  (ft) 

Layer  3 
-11.4 

Depth  to  bottom  (ft) 

-19.4 

cL( ft/sec) 

9000 

cu( ft/sec) 

9000 

cv(ft/sec) 

9000 

Poisson's  ratio 

0.25 

t0(bar£) 

51 

tan  (J> 

0.75 

YMAX  {kbars> 

2.1 

Tension  limit  (bars) 

-68 

MCP-02  with  the  Prandtl -Reuss  flow  rule  only 
MCP-02  with  no  yield  allowed 

MCP-02  with  yield  surface  shifted  only  before  yield  calculate 
MCP-02  with  yield  surface  shifted  only  after  yield  calculation 

MCP-09 : MC  2.12  with  the  changes 

Layer_i  Layer  2 Layer  3 Layer  4 
cu( ft/sec)  1250  11000  12000  11000 

cv(ft/sec)  500  4775  6164  4775 


101 


AFWL-TR-75-88 


MCP-12:  MCP-02  with  the  changes 


Layer  3 


Depth  to  top  (ft) 

-16.2 

cL (ft/ sec) 

6000 

cu(ft/sec) 

6000 

cv(ft/sec) 

6000 

Poisson's  ratio 

0.45 

tQ  (bars) 

34 

tan  <p 

0.75 

vmx  (bars) 

34 

Tension  limit  (bars) 

0 

MCP-21 : MCP-02  with  the  change  that  in  Layer  2 was  0.21  kbars 
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APPENDIX  B 
THE  AFT0N-2A  CODE 

The  theory  of  the  AFT0N-2A  code  has  been  documented  in  Air  Force  Weapons 
Laboratory  technical  reports  (Trulio,  1966;  Trulio,  et  a 1 . , 1969;  Niles,  et  al . , 

1971).  Since  the  work  accomplished  during  this  study  did  not  include  any  modi- 
fication of  that  theory,  this  description  is  not  significantly  different  from 
the  descriptions  provided  in  those  reports.  However,  in  order  to  provide  a 
more  complete  document,  the  description  of  the  theory  of  the  AFT0N-2A  code  by 
Niles,  et  al.,  1971  will  be  repeated. 

1.  FINITE  DIFFERENCE  MESHES  AND  ZONES  IN  AFT0N-2A 

The  finite  difference  technique  used  in  the  AFTON  codes  is  of  the  "time- 
marching" kind;  that  is,  the  space  continuum  is  replaced  by  a discrete  mesh  of 
points.  Starting  with  a system  in  a known  state  at  some  initial  time,  the 
variables  of  the  motion  are  updated  by  a discrete  time  increment  at  all  points 
of  the  space  mesh,  according  to  the  finite  difference  equations  of  motion.  The 
updating  process  is  then  repeated  using  the  just-calculated  values  of  the 
variables  of  the  motion  as  fresh  initial  value  data.  Owing  to  the  assumed 
symmetry  of  the  motion,  a space  mesh  for  AFT0N-2A  need  only  be  defined  as  an 
array  of  points  in  a single  azimuthal  plane,  the  variables  of  the  motion  having 
identical  values  at  corresponding  points  of  all  azimuthal  planes.  The  points  of 
an  AFT0N-2A  finite  difference  mesh  are  logically  equivalent  to  the  corner  points 
of  a set  of  unit  squares  which  cover  a rectangular  region  in  one-to-one  fashion. 

The  mesh  points  are  therefore  the  vertices  of  quadrilaterals  which  can  be  pro- 
duced by  the  continuous  distortion  of  a rectangular  array  of  unit  squares.  The 
region  of  two-dimensional  axisymmetric  flow  is  thus  covered  by  elementary 
quadrilaterals;  these  quadrilaterals  are  the  "zones"  of  the  finite  difference 
mesh.  Actually,  it  is  basic  to  the  method  of  differencing  which  underlies  the 
AFTON  codes  that  real  physical  systems  have  finite  extension  in  a direction 
normal  to  the  symmetry  plane  in  which  the  quadrilaterals  lie.  Thus,  a quadri- 
lateral zone  is  just  a cross  section  of  a quadrilateral  wedge  in  a single  azi- 
muthal plane  of  symmetry.  The  quadrilateral  wedge,  a solid  figure,  is  the 
basic  geometric  entity  of  the  AFT0N-2A  finite  difference  mesh  and  is  shown 
schematically  in  figure  B-l . It  is  a polyhedron  bounded  by  two  nearly  parallel 
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Figure  B-l.  Quadrilateral  Wedge 


azimuthal  planes  with  similar  quadrilateral  cross  sections  in  all  azimuthal 
planes  between  the  two.  Four  trapezoidal  faces  normal  to  the  central  quadri- 
lateral complete  the  polyhedron. 

The  integral  equations  and  associated  finite  difference  equations  which  are 
basic  to  AFT0N-2A  have  been  written  in  sufficient  generality  to  include  noi- 
Lagrangian  as  well  as  Lagrangian  descriptions  of  continuum  motion.  The  code 
itself  contains  a subroutine  which  defines  the  coordinate  system  to  be  used  for 
any  given  problem.  However,  the  Lagrangian  case  will  be  discussed  because  the 
finite  difference  technique  as  it  applies  to  AFT0N-2A  is  most  simply  explained 
for  that  case.  The  points  of  the  finite  difference  mesh  are  then  mass  points 
whose  velocities  provide  a discrete  approximation  to  the  material  velocity 
field  of  the  continuous  medium. 
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In  the  Lagrangian  case  each  quadrilateral  wedge  is  a finite  mass  element 
consisting  of  the  same  material  particles  at  one  time  as  at  any  other  time, 
and  a quadri later ial  zone--a  cross  section  of  a quadrilateral  wedge  in  a 
symmetry  plane--is  defined  by  one  specific  set  of  coplanar  particles.  Motion 
of  the  vertices  of  a quadrilateral  wedge  therefore  produces  a distortion  or 

strain  in  the  wedge  and  causes  changes  in  all  the  flow  variables  for  a finite 
element  of  material. 

2.  THE  CALCULATION  OF  THERMODYNAMIC  VARIABLES  IN  AFT0N-2A  FOR  LAGRANGIAN  MESHES 

The  variables  of  the  motion  are  divided  into  two  classes;  namely,  those 
associated  with  the  vertices  of  zones,  and  those  associated  with  their 
interiors.  The  first  class  (dynamic  variables)  consists  of  mesh  or  grid  point 
positions  and  their  time  derivatives  (i.e.,  their  velocities)  while  the  second 
class  (thermodynamic  variables)  includes  strain,  stress,  and  internal  energy. 

For  the  calculation  of  zone-centered  variables  two  assumptions  are  made. 

a.  A material  element  which  initially  occupies  the  region  enclosed  by  a 
quadrilateral  wedge  always  has  the  shape  of  a quadrilateral  wedge,  i.e., 
straight  lines  of  mass  points  deform  to  straight  lines  of  mass  points. 

b.  Zone-centered  variables  are  constant  in  value  throughout  a quadrilateral 
wedge  at  any  given  time,  and  also  change  at  a constant  rate  during  any 
particular  timestep. 

With  respect  to  assumption  (a),  the  particles  initially  comprising  a side 
of  a quadrilateral  zone  will,  in  general,  not  remain  col  inear;  likewise,  the 
corresponding  face  of  the  quadrilateral  wedge  associated  with  the  zone  usually 
will  not,  in  physical  reality,  remain  a quadrilateral.  Rather,  a trapezoidal 
Lagrangian  surface  of  the  quadrilateral  wedge  will  deform  into  a more  general 
curved  shape.  Assumption  (a),  therefore,  imposes  a nonphysical  constraint  on 
the  system,  which  is  part  of  the  price  paid  for  replacing  the  space  continuum 
by  a discrete  mesh  of  points,  Obviously,  assumption  (b)  entails  a similar 
nonphysical  restriction;  real  physical  stresses  and  strains  generally  vary 
continuously  over  finite  distances. 

The  calculation  of  the  change  in  the  volume  of  a quadrilateral  wedge  pro- 
duced by  the  motion  of  the  vertices  of  its  associated  quadrilateral  zone  pro- 
vides the  key  to  the  construction  of  the  finite  difference  equation  of  AFT0N-2A. 

In  making  the  calculation,  the  following  definitions  and  conventions  are 
adopted. 
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a.  V,  R,  U,  A denote  volume,  position  vector,  material  velocity  vector, 
and  vector  area,  respectively. 

b.  The  superscripts  n and  n-1  refer  to  a "later  time"  tn,  and  an  "earlier 
time"  tn’1,  separated  by  the  interval  4 tn'^^  = tn  - tn"^. 

c.  If  no  superscript  is  attached  to  a variable,  it  is  understood  to  be 
defined  at  some  time  between  tn~^  and  tn.  In  particular,  the  position 
vector  of  a point  without  a superscript  is  by  definition  equal  to  the 
arithmetic  mean  of  the  positions  of  the  point  at  the  two  times  tn  and 
t"'1,  i.e., 

r - 1/2  (rn  + rn_1  ) (6-1) 

n I 

d.  The  particle  velocity  of  a point  is  related  to  its  position  rn"  and 
j^n  at  the  times  tn"^  and  tn  according  to 

„ rn  - r"'1 

“ 4t"‘1/2  (B-2) 


e.  Position  and  velocity  subscripts  refer  to  the  mesh  points  labeled  as 
numbers  in  figures  B-l  and  B-2. 

The  underlined  subscripts  2,  3,  a^,  and  d^,  shown  schematically  in  figure  B-2, 
refer  to  points  on  the  side  of  zone  "a".  The  coordinates  of  point  j!,  for  ex- 
ample, are  defined  by  the  equations 


x2 


/ 2 2 \i/2 

| x2  + x2  x3  + x3 | 


(B-3) 


and 


v,  = ,(X3.  - xj  y?  + (x^.  - x2jy3 
x3  - x2 


(B-4) 


The  coordinates  of  the  points  3,  a_  and  d_  are  found  in  a similar  fashion. 

Equation  (B-2)  involves  the  kind  of  discretization  error  entailed  in  assump- 
tion (b)  above;  in  this  case,  the  velocity  is  taken  to  be  constant  over  a finite 
time  interval;  namely,  At.  An  exact  calculation  of  the  volume  of  a quadri- 
lateral wedge  shows  that 
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WMi 


1 


Vn  - Vn’^ 


At 


n-1/2 


4 

-SC 

i = l 


1/2 


u-1/2 


(B-5) 


where  the  index  i refers  to  the  vertices  1,2, 3, 4 of  the  quadrilateral  designated 
"a"  in  figure  B-2. 

In  the  limit  of  an  infinitesimal  timestep,  the  vector  area  is  composed  of 
a portion  of  the  two  trapezoidal  surfaces  whose  intersection  contains  the 
point  i.  For  example,  A3  corresponds  to  the  shaded  area  in  figure  B-l  when 

i = 3 and  thus 


A.3  = A_2_3  + A3_3_  (B  6) 

and  Asi  are  calculated  as  described  below  for  the  general  vector  area. 


The  vector  area  A . . is  the  trapezoidal  surface  of  the  quadrilateral  wedge 
between  the  vertices  i,  j of  figure  B-l.  The  sense  of  the  vector  area  Afj 
is  that  of  the  outer  normal  to  the  surface.  Thus,  for  example,  if  one  encoun- 

. i _ • 4-  ^ n 1 1 a rl  V'  n - 


lb  Lllau  Ul  UIIC  I 1W  — 

ters  point  i = 3,  and  then  the  point  j = 4,  as  the  perimeter  of  the  quadri- 
lateral  wedge  is  traversed  clockwise,  then 


A,j  - A„  » 1/2  [(13  - £*)  « (lo,  - 1,)] 

which  per  angle  <f>  can  be  shown  to  reduce  to 


(B-7 ) 


A34  = 1/2  ( x 3 + x4) 


y 3 - y 4 

X 4 - x3 
0 


(B-8) 


However,  the  vector  areas  A.  can  be  expressed  directly  in  terms  of  the  vertices 
of  quadrilateral  "a"  as  shown  below  for  A3 


A3X 

y?.(x2  + x3) 

- y4(X3  + X4)  + A2  3 + A 3 4 

a3  = 

A3y 

= 1/6 

-x2(x2  + x3) 

+ X4(X3  + x4) 

.0  . 

0 

- 

(B-9) 


k 


108 


__  - iiiiiMKH  . 


AFWL-TR- 75-88 


where 


A 2 3 = y"  - y?  (B-10) 

Equation  (B-5)  has  the  geometric  interpretation  that  the  change  in  the 
volume  of  a quadrilateral  wedge  in  a time  interval  At  is  equal  to  the  algebraic 
sum  of  the  volumes  swept  out  by  the  four  trapezoidal  faces  of  the  wedge  normal 
to  the  x-y  plane,  if  appropriate  portions  of  each  face  move  with  a uniform 
velocity  equal  to  the  velocity  of  the  vertex.  The  volume  change  so  calculated 
is  exact,  regardless  of  the  time  interval  At  or  of  the  positions  of  the  vertices 
of  the  quadrilateral  zone  at  the  beginning  and  end  of  the  interval. 

According  to  assumption  (b),  thermodynamic  variables  such  as  stresses  and 
internal  energies  are  considered  to  be  properties  of  quadrilateral  wedges  as  a 
whole.  These  variables  are  updated  for  general  stresses  and  strains  by  an 
extension  of  a standard  numerical  hydrodynamic  procedure  in  which  a finite 
difference  analog  of  the  First  Law  is  satisfied  simultaneously  with  the  con- 
stitutive equation  for  a given  medium.  In  the  hydrodynamic  case,  the  change 
in  the  internal  energy  of  a quadrilateral  wedge  is  just  its  volume  change 
given  by  equation  (B-5),  multiplied  by  the  negative  of  the  arithmetic  mean  of 
the  pressures  in  the  wedge  at  the  times  tn  ^ and  tn. 

En  - En_1  = -pq(vn  - Vn_1)  (B-l 1 ) 

where  PQ  denotes  (P  + Q)11"1^2.  If  an  equation  of  state  is  used  to  eliminate 
the  new  pressure  (i.e.,  the  pressure  at  time  tn)  from  the  finite  difference 
analog  of  the  First  Law,  then  the  fact  that  equations  of  state  generally  in- 
volve the  internal  energy  renders  the  First  Law  analog  an  implicit  equation  for 
the  new  internal  energy. 


Here  g,  the  equation  of  state,  is  some  (known)  function  of  two  variables,  and 
P,  E,  p denote  the  pressure,  internal  energy  and  density  of  the  quadrilateral 
wedge,  respectively,  the  mass  being  constant  in  the  Lagrangian  case  under  dis- 
cussion. Also,  Q is  a generalization  of  the  artificial  viscosity  of  von  Neumann 
and  Richtmyer  (1950).  Q is  computed  explicitly  knowing  V,  while  Pn  and  E 
should  be  obtained  by  solving  equations  (B-ll)  and  (B-12)  simultaneously. 
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In  the  calculation,  it  is  worth  noting  that  if  the  pressure  in  the  quadri- 
lateral wedge  were  indeed  uniform  and  equal  to  its  mean  value  on  the  time  in- 
terval At,  then  the  calculation  of  the  change  in  the  internal  energy  of  the 
wedge  as  well  as  its  volume  change  would  be  exact. 

For  general  axisymmetric  two-dimensional  motion,  the  procedure  for  writing 
an  exact  finite  difference  analog  of  the  First  Law  is  not  so  obvious  as  for 
hydrodynamic  motion.  In  fact,  an  exact  analog  of  the  First  Law  can  be  written 
only  for  triangular  zones  and  not  for  more  general  polygons  such  as  quadri- 
laterals. 


In  obtaining  a finite  difference  analog  of  the  First  Law  for  general  stress 
and  strain,  the  change  in  the  volume  of  the  zone,  as  given  in  equation  (B-5), 
is  of  prime  importance.  Introducing  this  expression  for  the  volume  change  into 
equation  (B-ll)  leads  directly  to  a finite  difference  analog  of  the  First  Law 
which  can  be  used  for  any  stress,  hydrodynamic  or  otherwise,  and  which  is 
exact  in  the  hydrodynamic  case  under  assumptions  (a)  and  (b).  This  combination 
of  equations  (B-3)  and  (B-ll)  for  the  zone  “a"  is 


En  . En-1 


= At 


I 

i = I 


“i 


li 


(B-l 3) 


where  for  hydrodynamic  motion,  the  forces  £1 , £4  in  equation  ( B-l 3)  are 

given  by  the  equations  of  the  form 


£1  = 


(B-l 4 ) 


To  compute  the  change  in  internal  energy  for  general  stresses,  the  scalar 

hydrodynamic  stress  (PQ)  of  equation  (B-l 4}  is  replaced  by  the  stress  tensor 

0...-  Again,  in  accord  with  assumption  (b),  a.,  is  assumed  to  be  constant  during 
"I  J ' J 

a timestep  throughout  any  particular  quadrilateral  wedge.  The  definitions  of 

the  forces  £1,  . . . , £4  then  become 

Fi  = a.  . • ( A + A . ) 

lj  \-ai  »d  J (B-l 5) 

etc.,  where  the  multiplication  called  for  in  equation  (B-l 5)  is  that  of  a matrix 
with  a vector. 
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As  an  expression  for  the  change  of  internal  energy,  the  right-hand  member 
of  equation  (B-13)  presents  one  obvious  problem:  its  terms  are  all  defined 

only  on  the  surface  of  a material  element,  whereas  "internal"  energy  is,  in 
fact,  a quantity  associated  in  an  essential  way  with  the  interior  of  a material 
region.  To  transform  equation  (B-13)  so  that  it  involves  only  interior  areas 
of  a quadrilateral  wedge,  an  elementary  geometric  theorem  is  used.  This 
theorem,  which  is  a cornerstone  of  the  finite  difference  method  embodied  in  the 
AFTON  codes,  simply  states  that  the  sum  of  the  vector  areas  of  any  polyhedral 
surface  is  zero,  where  the  sense  of  the  vector  area  associated  with  each  face 
of  the  polyhedron  is  understood  to  be  that  of  the  outer  normal  to  the  enclosed 
volume.  The  meaning  of  the  theorem  can  be  exhibited  in  the  following  geometric 
way.  Viewed  from  any  aspect  at  a sufficiently  great  distance,  a polyhedron 
presents  a cross  section  which  is  at  one  and  the  same  time  the  projection  of  the 
front  side  of  the  polyhedron  on  a plane  normal  to  the  viewer's  line  of  sight, 
and  also  of  its  back  side.  The  area  of  the  cross  section  is  equal  in  magnitude 
to  the  component  of  the  resultant  vector  area  of  the  plane  surfaces  making  up 
the  front  side  of  the  polyhedron,  and  is  also  the  negative  of  the  corresponding 
component  of  the  resultant  area  of  the  faces  of  the  back  side  (see  figure  B-3). 
Since  the  faces  of  the  front  and  back  side  make  up  the  entire  (closed)  poly- 
hedral surface,  the  sum  of  all  the  vector  areas  is  plainly  zero. 

With  respect  to  the  calculation  of  internal  energy  changes,  equation  (B-13) 
can  be  transformed  so  that  its  forces  refer  only  to  trapezoidal  surfaces  in  the 
interior  of  the  quadrilateral  wedge  plus  the  wedge  faces.  The  theorem  just 
discussed  implies,  for  example,  that  the  sum  of  the  trapezoidal  area  k 

^da’  Aaa/  anc^  1 * ^us  the  sum  0 f the  areas  of  the  two  azimuthal  surfaces  of 
the  wedge,  A^,  is  zero.  Therefore, 


4., 


+ 4.d  = ‘("da  + *aa  * 2A\ 


(B-16) 


The  vector  area,  , of  the  quadrilateral  cross  section  or  wedge  face  is 
defined  as  follows 


A =1/2 

— wi  ' 


- I.) 


la 


* ^a 


(B-17) 


in 
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Thus,  equation  (B-13)  can  be  written  in  the  form 


En  - En_1  = At 


faa(— 1 + ^ - Ui) 


da' 


+ F2a  ( U.3  - U2)  + FS3(U,  - U3) 


3a 


+ fwl  • Ui  + fw2  • u2  + ^3  • u3  + fW4  • U, 


(B-18) 


Equation  (B-18)  can  be  interpreted  as  a sum  of  internal  energy  changes  produced 
by  the  extension  of  material  in  directions  normal  to  the  forces  exerted  on 
specific  interior  or  wedge  face  surfaces  of  the  quadrilateral  wedge. 

3.  CALCULATION  OF  MOMENTUM 

In  addition  to  updating  the  thermodynamic  variables  of  a zone,  new  veloci- 
ties and  positions  of  the  mesh  points  also  need  to  be  calculated.  The  procedure 
used  to  update  velocities  is  based  on  the  principle  of  momentum  conservation, 
as  applied  to  a spatial  region  known  as  a "momentum  zone." 

As  in  the  case  of  velocities  and  positions,  momentum  zones  are  centered  at 
mesh  points.  The  momentum  zone  associated  with  a mesh  point  is  comprised  of  a 
precisely  defined  portion  of  each  of  the  four  thermodynamic  zones  which  share 
the  mesh  point  as  a common  vertex.  A thermodynamic  zone  is  therefore  divided 
into  four  pieces  each  of  which  is  associated  with  one,  and  only  one,  vertex 
for  the  purpose  of  the  momentum  calculation.  The  division  is  made  by  joining 
an  interior  point  of  the  zone,  called  its  "mid-point,"  to  certain  points  of  its 
edges;  for  example,  the  point  labeled  a in  figure  B-l  is  connected  to  the 
points  a_,  d,  2 and  3_.  The  trapezoidal  surfaces  bounded  by  pairs  of  points  such 
as  (a,aj  or  (a ,2)  (shaded  in  figure  B-4)  represent  a major  portion  of  the 
interior  areas  upon  which  the  stresses  are  imposed.  The  remaining  portion  of 
the  area  acted  upon  by  stress  is  subtended  by  the  wedge  faces  of  the  zone.  The 
momentum  zone  contains  a mass  of  material  equal  to  the  sum  of  a precise  portion 
of  each  of  the  four  thermodynamic  zones  which  have  as  a common  vertex  the  point 
about  which  the  momentum  zone  is  centered  (i.e.,  the  shaded  area  in  figure  B-5). 
The  momentum  zone  in  AFT0N-2A  is  then  a polyhedron  of  ten  faces.  Forces  exerted 
on  the  eight  trapezoidal  faces  and  two  wedge  faces  of  the  momentum  zone  produce 
an  acceleration  of  the  momentum  mass.  If  assumption  (b)  is  true  for  both  the 
forces  acting  on  the  momentum  zone  and  the  velocity  of  the  zone  during  an  entire 
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timestep,  then  the  change  in  momentum  during  a timestep  may  be  calculated 


exactly.  The  momentum  Mn  at  time  tn  can  now  be  updated  from  its  value  at 


n-1 


time  t . For  the  momentum  zone  conservation  of  momentum  is  expressed  by  the 
equation 


Mn  - Mn_1  = At 


F + F , + F + F , 
— ia  —2b  — sc  — 4d 


(B-l 9) 


In  the  above  expression,  the  force  F , is  related  to  three  of  the  forces  which 

1 a 


appear  in  equation  (B-18),  as  follows 


F = F + F , + F 

—Id  “del  — dd  — W1 


or 


£,a  = 


On(Aix  - a)  + o 12  A + a33  a 


O i 2 


(A,x  - a) 


+ a2  2 A 


iy 


(B-20) 


and  similarly  for  F ,,  F and  F ,.  The  sense  of  the  forces  F . . F „ is 
illustrated  in  figure  B-6.  a is  the  scalar  area  of  the  wedge  face  whose 
vertices  are  1,  c[,  a,  a^  The  velocity  of  the  mesh  point  on  which  the  momentum 
zone  is  centered  is  related  to  the  momentum  by  the  equation 


U"  = 5- 


m 


(B-21) 


where  m is  the  mass  of  the  momentum  zone.  A forward  extrapolation  in  time  is 
used  to  advance  the  velocity  from  one  timestep  to  the  next. 


yn+l/2  _ 2un  _ yn-1/2 


= 2Mf  _ yn-1/2 


(B-22) 
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4.  CALCULATION  OF  STRAIN 

The  discussion  of  the  calculation  of  strain  will  be  initially  restricted  to 
considering  an  axially  symmetric  wedge  (axis  of  symmetry  is  the  y-axis)  whose 
cross  section  in  the  x-y  plane  is  a triangle.  The  calculation  of  strain  for 
wedges  with  quadrilateral  cross  sections  will  be  treated  later.  Consider  a 
triangular  zone  with  vertices  1,  2,  3 in  its  unstrained  configuration  which, 
under  the  influence  of  external  forces,  is  strained  to  a new  configuration 
(figure  B-7).  Then,  in  axial  symmetry,  the  linear  transformation  which  takes 
a point  (x,y,z)  in  the  unstrained  state  to  the  strained  state  x',  y',  z')  is 
given  by 

x'  = an  x + ai2y 

y'  = a2i  x + a22y 

z'  = a33  z 

In  general  for  any  point  (x,  y,  z),  define 

“i 5 (*  - *,■)  ai = (x' 

1 6i5(y' 

Now  form  four  equations  in  the  four  unknowns  an,  ai2,  a2i»  a22>  which  are  ele- 
ments of  the  point-transformation  matrix  A1,  namely, 

oij  = an  ai  + a i2  Bi 

Oj  = 311  “2  + a12  ^2 

$ J - a2 i a3  + a22  Bi 

3^  = a2i  a2  + a22  62  (B-24) 


(B-23) 


from  which 
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Figure  B-7.  Schematic  of  an  Unstrained  and  Strained 
Triangular  Zone 


, (al32  - a23i)  • a12  - (aioti  - a2a2) 

du  - 5 


where  A 0 = 011B2  - a23i- 
ap 

Using  the  elements  of  the  transformation  matrix  A,  a new  matrix,  T,  can  be 
formed 


'a2  + a2 
11  12 

a 1 1 a 2 1 + a 12a22 

1 1 1 

1 1 2 

T = 

a 1 ia2 1 + a 1 2822 

a2  + a2 

1 1 2 

t22 
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The  eigenvalues  Xi  and  A2  of  T are  now  related  to  the  principal  extensions 
Ei  and  E2,  as  follows: 


Ei 


(B-27) 


where 


(til  + tl 2 ) ± 


/(tl  1 - 1 2 2 )2  + 4tj2| 


(B-28) 


The  third  principal  extension  E3  is  found  from  the  ratio  of  the  strained  to 
unstrained  volumes  according  to  the  equation 


f-  = E i E 2 E 3 


(B-29) 


The  three  principal  strains  are  related  to  the  extensions  by 


e 


x 


= Ei  - 1 


(B-30) 


The  principal  strain  axes  are  found  from  the  eigenvectors  of  T,  namely,  Ai , 
A2,  where 


Ai  = [(tn  - A)2  + t22 

-/2 

CNl 

+-> 
■ 1 

tii  - Aj 

r „ i 

1/2 

r 1 

A2  =[(t22  - A)2  + t22 

t2  2 - A 

CM 

4-> 

1 

Of  the  four  possible  vectors  which  can  be  formed  when  Ai  and  X2  are  substituted 
in  the  above  equations,  the  two  with  the  largest  positive  x-  and  y-  components 
are  chosen  as  the  principal  strain  axes. 

For  a wedge  whose  cross  section  is  a quadrilateral,  the  strain  cannot  be 
defined  uniquely.  The  convention  used  in  AFTON  to  define  the  strain  employs  an 
averaging  technique.  Each  diagonal  of  the  quadrilateral  divides  it  into  two 
triangles.  Elements  of  A-matrices  are  found  for  each  of  the  four  possible 
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triangles  formed  by  the  two  diagonals.  Then,  the  elements  of  T for  the  entire 
quadrilateral  are  formed  by  averaging  the  values  of  the  elements  of  the  A-matrix 
obtained  for  each  of  the  four  triangles. 

5.  CALCULATION  OF  STRESS 

Knowing  the  principal  strain  axes,  defined  by  the  eigenvectors  A and  A , 

* y 

the  principal  strains  can  then  be  rotated  into  the  laboratory  coordinate  system. 
The  components  of  strain  in  the  laboratory  coordinate  system  are  then  given  by 
the  relations 


£ xxL 
exyL 
eyyL 
ezzL 


(B-32) 


where  the  subscript  L indicates  that  the  strains  are  Lagrangian,  i.e.,  they  are 
computed  at  time  n for  the  mass  elements  that  occupied  generalized  coordinate 
cells  at  time  n-1.  Total  strain  increments  are  then  formed  as  follows: 


e 


n-i 

ij 


The  deviatoric  strain  increment  is  calculated  from  the  equation 


(B-33) 


A e 


1J 


= A£ij 


- 1/3  5 


1J 


Aeii 


( B- 34 ) 


Next,  the  Lagrangian  compression  n11  is  calculated  from  the  known  mass  and 
volume  of  this  material  element,  according  to  the  equation 


n 


n 

L 


m 


n-i 


V 


n 

L 


(B-35) 


The  compression  increment  AnL  = n"  - nj1"1.  and  the  excess  compression 
UL  = nL  ' ^ are  t*ien  f°rmed-  ^rom  the  equation  of  state  for  the  Lagrangian  mass 
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element,  the  new  mean  stress,  Pn,  its  derivative  K = dP/du,  and  the  shear 
modulus,  G,  are  computed.  The  deviatoric  elastic  stress  tensor  and  its  second 
invariant  are  then  computed  as  follows: 


(°t. j)"  ' Ho)""  - 2GAcij 


ne 1 _ , /9  e1  e1 
J2  - 1/2  0)j  o(j 


and  a function  Y is  evaluated  accora'ing  to  the  expression 

( ai  + oi2  P + a3  F ) 

Y = min 


where  o^.  and  k are  constants  and  P is  the  average  mean  stress  (Pn  + Pn-1)/2; 
the  yield  surface  equation  is  jf'  = Y2.  If  jf'<  Y2 , then 


rr'n  - „e' 

o • ■ - o . . 

U (E 

However,  if  J 2 > Y2,  then  the  incremental  mean  and  deviatoric  stresses  are 
formed  according  to  the  incremental  plastic  stress-strain  equations; 


KAu-  K(;£Ki  -fe 
' * (f  f I 


a!  - (a!.  Ac!.)  / dY \ . D 

j!  = G 2Ae ! • - — ^ Ul  _ a!  p 

J 2 ^2  Y — 


The  deviatoric  stress  is  then  computed  as  follows: 


_ I f I _ I N- 1 . . , 

°ij  - °ij  + 4aij 
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The  calculation  of  the  updated  components  of  Lagrangian  stress  is  completed 
using  the  equation 


6. . P' 
ij 


+ o! 


1J 


(B-43) 
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AUXILIARY  COMPUTATION  ROUTINES 


In  addition  to  the  AFTON-2A  code,  three  computer  programs  were  important 
to  the  results  and  displays  contained  in  this  document.  The  first  was  the 
ballistic  extension  model  used  to  extend  the  calculated  velocity  conditions  at 

16.4  msec  in  the  MIXED  COMPANY  II  simulation  to  the  final  grid  position  at 

616.4  msec.  The  second  was  a flow-field- parameter  display  program  that  used 
the  restart  dump  tapes  to  construct  spatial  displays  at  a particular  simulated 
time.  The  final  program  was  the  time-history  display  program  which  used  the 
target-point  data  tapes  to  construct  the  time  histories.  This  third  program 
involved  only  data  handling  and  simple  manipulations  and  will  not  be  described. 

However,  the  first  two  programs  included  calculational  procedures  that  are  im- 
portant to  their  use. 

1.  BALLISTIC  EXTENSION  ROUTINE 

The  ballistic  extension  model  had  to  describe  what  grid  points  were  moving, 
how  that  motion  changed  position  and  velocity  values  during  one  time  increment, 
and  what  grid  points  stopped  moving  during  the  time  increment.  The  major 
assumptions  of  the  routine  were  (1)  that  all  motion  occurred  under  the  influence 
of  gravity  only,  and  (2)  that  only  large  displacements  were  of  interest. 

A parameter,  S^,  was  used  at  each  grid  point  to  differentiate  between 
moving  and  stopped  grid  points.  For  moving  grid  points  the  value  of  this 
parameter  was  zero,  and  for  stationary  grid  points,  St  was  set  to  one.  At  the 
beginning  of  the  ballistic  extension  to  MC  2.12  the  grid  points  above  20  feet 
depth  and  within  a 35-foot  range  from  the  vertical  axis  were  considered  moving, 
while  all  other  grid  points  were  considered  stopped.  This  region  was  defined 
because  of  the  primary  interest  in  the  crater  region  and  because  the  velocity 
conditions  outside  that  region,  except  in  the  soil  layer,  would  result  in 
calculated  displacements  of  less  than  1 inch  under  ballistic  conditions. 

The  motion  of  each  grid  point  inside  that  region  was  calculated  at  time 
increments,  At,  of  2 milliseconds  until  three  conditions  were  simultaneously  > 

satisfied.  The  first  condition  was  that  at  least  one  of  the  grid  points  which 
were  originally  either  immediately  below  or  radially  away  from  the  grid  point 
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being  considered  had  stopped.  This  condition  was  met  if  the  value  of  at 
either  of  those  neighboring  grid  points  was  one.  The  second  condition  was  that 
the  vertical  velocity  component,  calculated  during  the  previous  time  increment, 
was  not  positive.  The  third  condition  was  that  the  material  density  in  the 
lower,  outward  quarter  volume  associated  with  the  grid  point  was  at  least 
1.5  gm/cc.  This  density  was  determined  by  dividing  the  mass  in  that  quarter 
zone,  which  was  constant  because  of  Lagrangian  motion,  by  the  volume  of  the 
quarter  zone  calculated,  using  the  volumetric  subroutine  of  AFT0N-2A,  at  the 
beginning  of  the  time  increment.  If  all  three  of  these  conditions  were  met, 
the  value  of  at  the  zone  being  considered  was  set  to  one. 


The  motion  of  each  grid  point  was  then  calculated  based  on  the  value  of 

Sj..  If  the  value  of  at  the  grid  point  was  one,  both  the  horizontal  and 

vertical  velocity  components  were  set  to  zero  and  the  position  coordinates  of 

the  grid  point  remained  constant.  If  the  value  of  at  the  grid  point  was 

zero,  then  the  horizontal  velocity  component,  U , remained  constant  and  the 

T X 

vertical  velocity  component  U , at  the  end  of  the  time  increment  was  determined 
by 


Uy  = U°  - gAt  (C-l) 

o 

where  is  the  vertical  velocity  component  at  the  beginning  of  the  time 
increment  and  g is  the  gravitational  acceleration.  The  position  coordinates 

i i 

of  the  grid  point  (X  , Y ) at  the  end  of  the  time  increment  were  then 


X'  = X + U At 

A 


Y ' = Y + At 


0 Q ' 

Uy  - f it 


(C-2) 


where  (X  , Y ) are  the  position  coordinates  at  the  beginning  of  the  time 
increment. 


2.  FLOW  FIELD  DISPLAY  ROUTINES 


The  calculation  space  displays  were  used  to  display  calculated  motion  and 
thermodynamic  parameters  at  appropriate  positions.  These  values,  except  for 
the  maximum  pressure  contour  plot,  were  for  a particular  moment  of  simulated 
time  and  were  represented  by  vector  arrows  that  began  at  the  calculational 
position,  indicated  the  vector  direction,  and  had  lengths  scaled  to  the 
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magnitude  of  the  value.  A vector  would  not  be  drawn,  however,  if  (1)  the 
calculaf.ional  position  was  outside  the  display  field,  (2)  the  display  length 
of  the  vector  was  less  than  a minimum  value,  or  (3)  the  length  of  the  vector 
would  cause  it  to  extend  beyond  the  display  borders.  Minimum  length  vectors 
were  0.01  inch  on  the  velocity  vector  plots  and  0.03  inch  on  the  acceleration, 
principal  stress,  and  calculation  grid  plots.  The  displays  that  required 
significant  data  manipulation  were  the  plots  of  acceleration  vectors,  principal 
stress  axes,  and  maximum  pressure  contours. 

The  acceleration  vectors  and  principal  stress  vectors  were  represented  by 
arrows  with  lengths  proportional  to  the  square  root  of  magnitude  but  were 
constructed  to  maintain  the  exact  vector  direction.  The  second  condition 
required  that 


(C-3) 


where  1 and  1 are  the  x and  y components  of  the  display  vector,  and  Ax  and 
x y 

Ay  are  the  linearly  scaled  components  of  the  quantity  to  be  represented.  The 
first  condition  is  satisfied  if 

’5+  ’HAx  + Ay),/!  (C-4) 

Solving  these  two  relations  for  1 implied 


where  the  sign  for  1 was  chosen  to  be  the  same  as  A . The  value  of  1 was 
then  determined  from  relation  (C-3).  However,  this  scheme  would  not  work  on  a 
computer  if  A^  was  near  zero.  Since  the  Calcomp  hardware  will  not  plot  dis- 
tances less  than  0.01  inch,  an  alternate  scheme  was  used  when  Ay  was  less  than 
0.0001.  The  alternate  scheme  set  1 equal  zero  and  set 


1 


x 


+ 


(C-6) 
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where  the  sign  of  1 was  chosen  to  be  the  same  as  A . 

A A 

The  acceleration  vector  was  determined  from  the  velocity  and  timestep  infor- 
mation in  the  restart  dumps.  The  AFT0N-2A  code  retains  velocity  components 
(Unv,  Ll[| ) at  a simulated  time  and  velocity  components  (ll^+1//%  u"+  ^jextrapo- 
lated  one-half  timestep  ahead.  The  acceleration  components  |ax,  ay ) were 
determined  by 


ax  - 2 (u"/V>  . u"x)«t" 
ay  - 2 (u  J+'/2  - Uj)/At" 


where  Atn  is  the  timestep. 

The  calculation  of  principal  stresses  and  principal  stress  direction  was 
based  on  the  Mohr  circle  construction.  The  maximum  and  minimum  principal 
stresses  were  found  by 


CT.  . 0.5  (oxx  + oyy)  + R 

<”  ■ °'5  (°xx  + °yy)  - R (c-8> 

where  R = [o.25(axx  - ayy) 2 + and  axx,  ayy,  and  oyy  represent  the 

radial,  vertical,  and  shear  stresses  from  the  restart  dump.  The  maximum  princi- 
pal stress  was  considered  to  be  in  the  x principal  axis  direction  unless  oyy 
was  greater  than  a . The  angle  of  rotation,  0,  of  the  principal  axes  from  the 

A A 

page  coordinates  was 


0 = 0.5  tan-1 


2a 


xy_ 


laxx  " CTyyl 


(C-9) 


unless  | a | > 100  • |a  - a | in  which  case  0 = ±45°  with  the  sign  chosen 
\y  XX  yy 

the  same  as  a . A positive  0 would  rotate  the  principal  stress  axes  from  the 
page  axes  in  a clockwise  manner.  This  scheme  was  used  unless  the  stress  values 
indicated  the  material  had  separated  at  the  thermodynamic  point  being  consid- 
ered, in  which  case  an  X was  centered  at  the  point. 
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A maximum  pressure  contour  plot  involved  first  locating  the  original  posi- 
tion of  the  thermodynamic  grid  point  and  then  determining  the  maximum  pressure 
experienced  by  that  point.  The  original  position  of  the  grid  point  was  deter- 
mined by  first  subtracting  the  displacements  of  the  four  surrounding  mesh 
points  from  their  positions  at  the  time  of  the  restart  dump  to  determine  the 
original  grid  positions.  Then  the  excess  compression  value  was  equated  to  the 
maximum  excess  compression  value  of  the  thermodynamic  point  being  considered 
and  the  AFTON-^A  equatior.-of-state  routine  was  used  to  determine  the  maximum 
pressure  value.  These  values  were  then  contoured  with  the  use  of  a contour 
mapping  subroutine  which  was  provided  by  the  computer  support  division  of  the 
Air  Force  Cambridge  Research  Laboratories . 
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ATTN:  D.  H.  Scott 
ATTN:  E.  M.  Shoemaker 
ATTN:  L.  A.  Soderblom 
ATTN:  G.  A.  Swann 
3 cy  ATTN:  D.  J.  Roddy 

Department  of  the  Interior 
U.S.  Geological  Survey 

ATTN:  Daniel  J.  Milton 
ATTN:  Richard  S.  Pike,  Jr. 

ATTN:  Don  E.  Wilhelms 
ATTN:  Howard  G.  Wilshirc 

ATTN:  Cecil  B.  Raleigh,  Earthquake  Res.  Center 
ATTN:  John  H.  Healy 


OTHER  GOVERNMENT  AGENCIES  (Continued) 

Department  of  the  interior 
Bureau  of  Mines 

ATTN:  Dr.  Leonard  A.  Obert 

Bureau  of  Mines 
Twin  Cities  Research  Center 
ATTN:  Dr.  T.  C.  Atchison 
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U.  S.  Geological  Survey 

ATTN:  Edward  C.  T.  Chao 

DEPARTMENT  OP  DEFENSE  CONTRACTORS 

Aerospace  Corporation 

ATTN:  Dr.  Prc.n  Mather 

Aerospace  Corporation 

ATTN:  Mr.  Warren  Pfefferle 
ATTN:  l)r.  Mason  B.  Watson 

Agbabian  Associates 

ATTN:  Dr.  Mike  Agbabian 

Applied  Theory,  Inc. 

ATTN:  Dr.  J.  Trulio 
ATTN:  Neil  Pearl 

The  Boeing  Company 

ATTN:  Mr.  Ron  Carlson 

Brown  Engineering  Company,  Inc. 

ATTN:  Manu  Pate! 

California  Institute  of  Technology 
ATTN:  Dr.  Thomas  J.  Ahrens 
ATTN:  Dr.  Leon  T.  Silver 

California  Researeh  & Technology,  Inc. 

ATTN:  M.  Rosenblatt 

Civil/Nuclear  Systems  Corporation 
ATTN-  Dr.  Robert  Crawford 

1IT  Research  Institute 

ATTN:  Tech  Library 

Institute  for  Defense  Analyses 

ATTN:  Tech  information  Offiee 

institute  of  Geophysics  6 Planetary  Physics 
ATTN:  Orson  J.  Anderson 

General  Electric  Company 

TEMPO-Center  for  Advanced  Studies 
ATTN:  DASLAC 

Lockheed  Missiles  & Space  Company 

ATTN:  Dr.  Ronald  E.  Meyerott,  Dept.  50-01 
Bldg.  201 

The  Lunar  Spaee  Institute 

ATTN:  Dr.  Frederick  Horz 

Massachusetts  Institute  of  Technology 
ATTN:  Professor  Eugene  Simmons 
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DEPARTMENT  OF  DEFENSE  CONTRACTORS  (Continued) 


McDonnell- Douglas  Corporation 
ATTN:  Mr.  Ken  MeClymonds 

Merritt  Cases,  Inc. 

ATTN:  J.  L.  Merritt 

Occidental  College 

Department  ol'  Geology 

ATTN:  David  Cummings 

Paeifiea  Technology 

ATTN:  Dr.  R.  T.  Allen 
ATTN:  Dr.  R.  1,.  Bjork 

Physics  International  Company 

ATTN:  Doc.  Con.  lor  Dr.  Charles  Godfrey 
ATTN:  Doc.  Con.  for  Mr.  Fred  M.  Sauer 
ATTN:  Doc.  Con.  for  Dr.  Robert  Swift 
ATTN:  Doc.  Con.  for  Mr.  Dennis  Orphal 

R & D Associates 

ATTN:  Dr.  Harold  L.  Brode 
ATTN:  Mr.  Bob  Thompson 
ATTN:  Mr.  John  Levesque 
ATTN;  Dr.  C.  I>.  Knowles 

d cy  ATTN:  Dr.  Henry  Cooper 

Rand  Corporation 

ATTN:  Dr.  C.  C.  Mow 

Science  Applications,  hie. 

ATTN:  Mr.  Mike  McKay 

Science  Applications,  Inc. 

ATTN:  Bill  Layson 

Science  Applications,  Inc. 

ATTN:  Dr.  D.  Maxwell 
ATTN:  Mr.  R.  Hoffmann 

Stanford  Research  Institute 
ATTN:  Dr.  Carl  Peterson 
ATTN:  Mr.  George  Abrahamson 


Systems,  Science  £ Software 
ATTN:  Dr.  Ted  Cherry 
ATTN:  Dr.  Ronald  R.  Grine 
ATTN:  Dr.  D.  Riney 
ATTN:  Mr.  Bob  A.  Ken 
ATTN:  Document  Control 

Terra  Tek,  Inc. 

ATTN:  Dr.  II.  R.  Pratt 
ATTN:  Dr.  J.  N.  Johnson 

TRW  Systems  G roup 

ATTN:  Mr.  Bing  F'ay 

TRW  Systems  Group 

ATTN:  Mr.  Norm  Lipner 
ATTN:  Dr.  Peter  K.  Dai,  1U-2178 
ATTN:  Dr.  Benjamin  Sussholtz 
ATTN:  Mr.  Edgar  Wong 

University  of  Illinois 

ATTN:  Dr.  Nathan  M.  Newmark 
ATTN:  Dr.  Skip  Hcndron 
ATTN;  Dr.  Bill  Hall 

University  of  New  Mexico 

Civil  Engineering  Research  Facility 
ATTN:  Mr.  Del  Calhoun 
ATTN:  Mr.  D.  J.  Higgins 

University  of  Texas 

Department  of  Geological  Sciences 
ATTN:  William  R.  Muehlberger 

Weidlinger  Associates,  Consulting  Engineers 
ATTN:  Dr.  Melvin  L.  Baron 
ATTN:  Ivan  Nelson 
ATTN:  Ivan  Sandler 

Weidlinger  Associates,  Consulting  Engineers 
ATTN:  Dr.  J.  Isenberg 
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